{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CDEC: Unsupervised Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\admin-karim\\appdata\\local\\programs\\python\\python35\\lib\\site-packages\\matplotlib\\__init__.py:1405: UserWarning: \n",
      "This call to matplotlib.use() has no effect because the backend has already\n",
      "been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n",
      "c:\\users\\admin-karim\\appdata\\local\\programs\\python\\python35\\lib\\site-packages\\sklearn\\utils\\linear_assignment_.py:21: DeprecationWarning: The linear_assignment_ module is deprecated in 0.21 and will be removed from 0.23. Use scipy.optimize.linear_sum_assignment instead.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "from time import time\n",
    "from keras.datasets import mnist\n",
    "import numpy as np\n",
    "np.random.seed(10)\n",
    "import numpy as np\n",
    "import keras.backend as K\n",
    "from keras.engine.topology import Layer, InputSpec\n",
    "from keras.layers import Dense, Input\n",
    "from keras.models import Model\n",
    "from keras.optimizers import SGD\n",
    "from keras import callbacks\n",
    "from keras.initializers import VarianceScaling\n",
    "from sklearn.cluster import KMeans\n",
    "import metrics\n",
    "\n",
    "from keras.models import Model\n",
    "from keras import backend as K\n",
    "from keras import layers\n",
    "from keras.layers import Input, Dense, Conv2D, MaxPooling2D, UpSampling2D, Flatten, Reshape, Conv2DTranspose\n",
    "from keras.models import Model\n",
    "import numpy as np\n",
    "import os \n",
    "from keras.preprocessing.image import load_img\n",
    "import _pickle as cPickle\n",
    "import _pickle \n",
    "import seaborn as sns\n",
    "import sklearn.metrics\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.manifold import TSNE\n",
    "%matplotlib inline\n",
    "import gzip\n",
    "import numpy as np\n",
    "from PIL import Image\n",
    "import matplotlib\n",
    "import os\n",
    "# For plotting graphs via ssh with no display\n",
    "# Ref: https://stackoverflow.com/questions/2801882/generating-a-png-with-matplotlib-when-display-is-undefined\n",
    "matplotlib.use('Agg')\n",
    "from keras.preprocessing.image import load_img\n",
    "from matplotlib import pyplot as plt\n",
    "from numpy import float32\n",
    "from sklearn import metrics\n",
    "from sklearn.cluster.k_means_ import KMeans\n",
    "from sklearn import manifold\n",
    "from sklearn.utils.linear_assignment_ import linear_assignment\n",
    "from sklearn import preprocessing\n",
    "from sklearn.utils import shuffle\n",
    "from keras.layers import BatchNormalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convolutional DEC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: convolutional auto encoder using `Conv2DTranspose`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autoencoderConv2D_1(input_shape=(512, 512, 1), filters=[32, 64, 128, 10]):\n",
    "    input_img = Input(shape=input_shape)\n",
    "    if input_shape[0] % 8 == 0:\n",
    "        pad3 = 'same'\n",
    "    else:\n",
    "        pad3 = 'valid'\n",
    "    x = Conv2D(filters[0], 5, strides=2, padding='same', activation='relu', name='conv1', input_shape=input_shape)(input_img)\n",
    "\n",
    "    x = Conv2D(filters[1], 5, strides=2, padding='same', activation='relu', name='conv2')(x)\n",
    "\n",
    "    x = Conv2D(filters[2], 3, strides=2, padding=pad3, activation='relu', name='conv3')(x)\n",
    "\n",
    "    x = Flatten()(x)\n",
    "    encoded = Dense(units=filters[3], name='embedding')(x)\n",
    "    x = Dense(units=filters[2]*int(input_shape[0]/8)*int(input_shape[0]/8), activation='relu')(encoded)\n",
    "\n",
    "    x = Reshape((int(input_shape[0]/8), int(input_shape[0]/8), filters[2]))(x)\n",
    "    x = Conv2DTranspose(filters[1], 3, strides=2, padding=pad3, activation='relu', name='deconv3')(x)\n",
    "\n",
    "    x = Conv2DTranspose(filters[0], 5, strides=2, padding='same', activation='relu', name='deconv2')(x)\n",
    "\n",
    "    decoded = Conv2DTranspose(input_shape[2], 5, strides=2, padding='same', name='deconv1')(x)\n",
    "    return Model(inputs=input_img, outputs=decoded, name='AE'), Model(inputs=input_img, outputs=encoded, name='encoder')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolutional auto encoder using `UpSampling2D`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = VarianceScaling(scale=1. / 3., mode='fan_in', distribution='uniform')\n",
    "\n",
    "def autoencoderConv2D_2(img_shape=(508, 508, 3)):\n",
    "    \"\"\"\n",
    "    Conv2D auto-encoder model.\n",
    "    Arguments:\n",
    "        img_shape: e.g. (512, 512, 1) for MNIST\n",
    "    return:\n",
    "        (autoencoder, encoder), Model of autoencoder and model of encoder\n",
    "    \"\"\"\n",
    "    input_img = Input(shape=img_shape)\n",
    "    # Encoder\n",
    "    x = Conv2D(128, (3, 3), activation='relu', padding='same', strides=(2, 2))(input_img)\n",
    "    x = Conv2D(64, (3, 3), activation='relu', padding='same', strides=(2, 2))(x)\n",
    "    x = Conv2D(32, (3, 3), activation='relu', padding='same', strides=(2, 2))(x)\n",
    "    x = BatchNormalization()(x)\n",
    "    shape_before_flattening = K.int_shape(x)\n",
    "    # at this point the representation is (4, 4, 8) i.e. 128-dimensional\n",
    "    x = Flatten()(x)\n",
    "    encoded = Dense(4, kernel_initializer=init, activation='relu', name='encoded')(x)\n",
    "    x = BatchNormalization()(encoded)\n",
    "\n",
    "    # Decoder\n",
    "    x = Dense(np.prod(shape_before_flattening[1:]), activation='relu', kernel_initializer=init)(encoded)\n",
    "    # Reshape into an image of the same shape as before our last `Flatten` layer\n",
    "    x = Reshape(shape_before_flattening[1:])(x)\n",
    "\n",
    "    x = Conv2D(32, (3, 3), activation='relu', padding='same')(x)\n",
    "    x = BatchNormalization()(x)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    x = Conv2D(64, (3, 3), activation='relu', padding='same')(x)\n",
    "    x = BatchNormalization()(x)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    x = Conv2D(128, (3, 3), activation='relu')(x)\n",
    "    x = BatchNormalization()(x)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    decoded = Conv2D(3, (3, 3), activation='sigmoid', padding='same')(x)\n",
    "\n",
    "    return Model(inputs=input_img, outputs=decoded, name='AE'), Model(inputs=input_img, outputs=encoded, name='encoder')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pick a convolutional autoencoder (`autoencoderConv2D_1` or `autoencoderConv2D_2`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder, encoder = autoencoderConv2D_2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_3 (InputLayer)         (None, 508, 508, 3)       0         \n",
      "_________________________________________________________________\n",
      "conv2d_15 (Conv2D)           (None, 254, 254, 128)     3584      \n",
      "_________________________________________________________________\n",
      "conv2d_16 (Conv2D)           (None, 127, 127, 64)      73792     \n",
      "_________________________________________________________________\n",
      "conv2d_17 (Conv2D)           (None, 64, 64, 32)        18464     \n",
      "_________________________________________________________________\n",
      "batch_normalization_2 (Batch (None, 64, 64, 32)        128       \n",
      "_________________________________________________________________\n",
      "flatten_3 (Flatten)          (None, 131072)            0         \n",
      "_________________________________________________________________\n",
      "encoded (Dense)              (None, 4)                 524292    \n",
      "_________________________________________________________________\n",
      "dense_3 (Dense)              (None, 131072)            655360    \n",
      "_________________________________________________________________\n",
      "reshape_3 (Reshape)          (None, 64, 64, 32)        0         \n",
      "_________________________________________________________________\n",
      "conv2d_18 (Conv2D)           (None, 64, 64, 32)        9248      \n",
      "_________________________________________________________________\n",
      "batch_normalization_4 (Batch (None, 64, 64, 32)        128       \n",
      "_________________________________________________________________\n",
      "up_sampling2d_7 (UpSampling2 (None, 128, 128, 32)      0         \n",
      "_________________________________________________________________\n",
      "conv2d_19 (Conv2D)           (None, 128, 128, 64)      18496     \n",
      "_________________________________________________________________\n",
      "batch_normalization_5 (Batch (None, 128, 128, 64)      256       \n",
      "_________________________________________________________________\n",
      "up_sampling2d_8 (UpSampling2 (None, 256, 256, 64)      0         \n",
      "_________________________________________________________________\n",
      "conv2d_20 (Conv2D)           (None, 254, 254, 128)     73856     \n",
      "_________________________________________________________________\n",
      "batch_normalization_6 (Batch (None, 254, 254, 128)     512       \n",
      "_________________________________________________________________\n",
      "up_sampling2d_9 (UpSampling2 (None, 508, 508, 128)     0         \n",
      "_________________________________________________________________\n",
      "conv2d_21 (Conv2D)           (None, 508, 508, 3)       3459      \n",
      "=================================================================\n",
      "Total params: 1,381,575\n",
      "Trainable params: 1,381,063\n",
      "Non-trainable params: 512\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "autoencoder.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load breast biohistology images for convolutional input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loadDataset():\n",
    "    data= []\n",
    "    labels =[]\n",
    "    root ='/home/rkarim/Training_data/'\n",
    "\n",
    "    for rootName,dirName,fileNames in os.walk(root):\n",
    "        if(not rootName == root):\n",
    "            for fileName in fileNames:\n",
    "                imgGray = load_img(rootName+'/'+fileName,color_mode='rgb')\n",
    "                if rootName.split('/')[1] == 'Benign':\n",
    "                    labels+=[0]\n",
    "                elif rootName.split('/')[1]== 'InSitu':\n",
    "                    labels+=[1]\n",
    "                elif rootName.split('/')[1]  == 'Invasive':\n",
    "                    labels+=[2]\n",
    "                else:\n",
    "                    labels+=[3]\n",
    "                  \n",
    "                transformed=transform.resize(np.array(imgGray),(508,508))\n",
    "                data += [transformed.reshape((transformed.shape[0],transformed.shape[1],3))]        \n",
    "          \n",
    "    data = np.stack(data)\n",
    "    labels = np.stack(labels)\n",
    "    #data,labels = shuffle(data,labels,random_state = 0)\n",
    "    \n",
    "    return [data,labels]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(266, 508, 508, 3)"
      ]
     },
     "execution_count": 263,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from keras.datasets import mnist\n",
    "\n",
    "x, y =  loadDataset()\n",
    "x.shape\n",
    "\n",
    "#x = x.reshape(x.shape + (1,))\n",
    "#x = np.divide(x, 255.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \n",
    "     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, \n",
    "     1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, \n",
    "     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, \n",
    "     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, \n",
    "     3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, \n",
    "     3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.array(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 267,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_clusters = len(np.unique(y))\n",
    "x.shape\n",
    "n_clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pretrain covolutional autoencoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {},
   "outputs": [],
   "source": [
    "pretrain_epochs = 1\n",
    "batch_size = 128\n",
    "save_dir = 'results/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.compile(optimizer='adam', loss='mse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/1\n",
      "266/266 [==============================] - 248s 932ms/step - loss: 0.0619\n"
     ]
    }
   ],
   "source": [
    "autoencoder.fit(x, x, batch_size=batch_size, epochs=pretrain_epochs)\n",
    "autoencoder.save_weights(save_dir+'/conv_ae_weights.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [],
   "source": [
    "autoencoder.load_weights(save_dir+'/conv_ae_weights.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build clustering model with convolutional autoencoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ClusteringLayer(Layer):\n",
    "    \"\"\"\n",
    "    Clustering layer converts input sample (feature) to soft label, i.e. a vector that represents the probability of the\n",
    "    sample belonging to each cluster. The probability is calculated with student's t-distribution.\n",
    "\n",
    "    # Example\n",
    "    ```\n",
    "        model.add(ClusteringLayer(n_clusters=10))\n",
    "    ```\n",
    "    # Arguments\n",
    "        n_clusters: number of clusters.\n",
    "        weights: list of Numpy array with shape `(n_clusters, n_features)` witch represents the initial cluster centers.\n",
    "        alpha: degrees of freedom parameter in Student's t-distribution. Default to 1.0.\n",
    "    # Input shape\n",
    "        2D tensor with shape: `(n_samples, n_features)`.\n",
    "    # Output shape\n",
    "        2D tensor with shape: `(n_samples, n_clusters)`.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, n_clusters, weights=None, alpha=1.0, **kwargs):\n",
    "        if 'input_shape' not in kwargs and 'input_dim' in kwargs:\n",
    "            kwargs['input_shape'] = (kwargs.pop('input_dim'),)\n",
    "        super(ClusteringLayer, self).__init__(**kwargs)\n",
    "        self.n_clusters = n_clusters\n",
    "        self.alpha = alpha\n",
    "        self.initial_weights = weights\n",
    "        self.input_spec = InputSpec(ndim=2)\n",
    "\n",
    "    def build(self, input_shape):\n",
    "        assert len(input_shape) == 2\n",
    "        input_dim = input_shape[1]\n",
    "        self.input_spec = InputSpec(dtype=K.floatx(), shape=(None, input_dim))\n",
    "        self.clusters = self.add_weight((self.n_clusters, input_dim), initializer='glorot_uniform', name='clusters')\n",
    "        if self.initial_weights is not None:\n",
    "            self.set_weights(self.initial_weights)\n",
    "            del self.initial_weights\n",
    "        self.built = True\n",
    "\n",
    "    def call(self, inputs, **kwargs):\n",
    "        \"\"\" student t-distribution, as same as used in t-SNE algorithm.\n",
    "         Measure the similarity between embedded point z_i and centroid µ_j.\n",
    "                 q_ij = 1/(1+dist(x_i, µ_j)^2), then normalize it.\n",
    "                 q_ij can be interpreted as the probability of assigning sample i to cluster j.\n",
    "                 (i.e., a soft assignment)\n",
    "        Arguments:\n",
    "            inputs: the variable containing data, shape=(n_samples, n_features)\n",
    "        Return:\n",
    "            q: student's t-distribution, or soft labels for each sample. shape=(n_samples, n_clusters)\n",
    "        \"\"\"\n",
    "        q = 1.0 / (1.0 + (K.sum(K.square(K.expand_dims(inputs, axis=1) - self.clusters), axis=2) / self.alpha))\n",
    "        q **= (self.alpha + 1.0) / 2.0\n",
    "        q = K.transpose(K.transpose(q) / K.sum(q, axis=1)) # Make sure each sample's 10 values add up to 1.\n",
    "        return q\n",
    "\n",
    "    def compute_output_shape(self, input_shape):\n",
    "        assert input_shape and len(input_shape) == 2\n",
    "        return input_shape[0], self.n_clusters\n",
    "\n",
    "    def get_config(self):\n",
    "        config = {'n_clusters': self.n_clusters}\n",
    "        base_config = super(ClusteringLayer, self).get_config()\n",
    "        return dict(list(base_config.items()) + list(config.items()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {},
   "outputs": [],
   "source": [
    "clustering_layer = ClusteringLayer(n_clusters, name='clustering')(encoder.output)\n",
    "model = Model(inputs=encoder.input, outputs=clustering_layer)\n",
    "model.compile(optimizer='adam', loss='kld')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: initialize cluster centers using k-means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=n_clusters, n_init=5)\n",
    "y_pred_kmeans = kmeans.fit_predict(encoder.predict(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2819548872180451"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metrics.accuracy_score(y, y_pred_kmeans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_last = np.copy(y_pred_kmeans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.get_layer(name='clustering').set_weights([kmeans.cluster_centers_])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: deep clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 299,
   "metadata": {},
   "outputs": [],
   "source": [
    "#computing an auxiliary target distribution\n",
    "def target_distribution(q):\n",
    "    weight = q ** 2 / q.sum(0)\n",
    "    return (weight.T / weight.sum(1)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 300,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = 0\n",
    "index = 0\n",
    "maxiter = 1000\n",
    "update_interval = 10\n",
    "index_array = np.arange(x.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "metadata": {},
   "outputs": [],
   "source": [
    "tol = 0.01 # tolerance threshold to stop training"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Start training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 0: acc = 0.28195, nmi = 0.04127, ari = 0.01177  ; loss= 0\n",
      "Iter 10: acc = 0.27820, nmi = 0.00000, ari = 0.00000  ; loss= 0.21101\n",
      "Iter 20: acc = 0.27820, nmi = 0.00000, ari = 0.00000  ; loss= 0.0\n",
      "delta_label  0.0 < tol  0.01\n",
      "Reached tolerance threshold. Stopping training.\n"
     ]
    }
   ],
   "source": [
    "for ite in range(int(maxiter)):\n",
    "    if ite % update_interval == 0:\n",
    "        q = model.predict(x, verbose=0)\n",
    "        p = target_distribution(q)  # update the auxiliary target distribution p\n",
    "\n",
    "        # evaluate the clustering performance\n",
    "        y_pred = q.argmax(1)\n",
    "        if y is not None:\n",
    "            acc = np.round(metrics.accuracy_score(y, y_pred), 5)\n",
    "            nmi = np.round(metrics.mutual_info_score(y, y_pred), 5)\n",
    "            ari = np.round(metrics.adjusted_rand_score(y, y_pred), 5)\n",
    "            loss = np.round(loss, 5)\n",
    "            print('Iter %d: acc = %.5f, nmi = %.5f, ari = %.5f' % (ite, acc, nmi, ari), ' ; loss=', loss)\n",
    "\n",
    "        # check stop criterion\n",
    "        delta_label = np.sum(y_pred != y_pred_last).astype(np.float32) / y_pred.shape[0]\n",
    "        y_pred_last = np.copy(y_pred)\n",
    "        if ite > 0 and delta_label < tol:\n",
    "            print('delta_label ', delta_label, '< tol ', tol)\n",
    "            print('Reached tolerance threshold. Stopping training.')\n",
    "            break\n",
    "    idx = index_array[index * batch_size: min((index+1) * batch_size, x.shape[0])]\n",
    "    loss = model.train_on_batch(x=x[idx], y=p[idx])\n",
    "    index = index + 1 if (index + 1) * batch_size <= x.shape[0] else 0\n",
    "\n",
    "model.save_weights(save_dir + '/conv_DEC_model_final.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the clustering model trained weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.load_weights(save_dir + '/conv_DEC_model_final.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 304,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acc = 0.28000, nmi = 0.00000, ari = 0.00000  ; loss= 0.0\n"
     ]
    }
   ],
   "source": [
    "#Eval.\n",
    "\n",
    "q = model.predict(x, verbose=0)\n",
    "p = target_distribution(q)  # update the auxiliary target distribution p\n",
    "\n",
    "# evaluate the clustering performance\n",
    "y_pred = q.argmax(1)\n",
    "if y is not None:\n",
    "    acc = np.round(metrics.accuracy_score(y, y_pred), 2)\n",
    "    nmi = np.round(metrics.mutual_info_score(y, y_pred), 2)\n",
    "    ari = np.round(metrics.adjusted_rand_score(y, y_pred), 2)\n",
    "    loss = np.round(loss, 5)\n",
    "    print('Acc = %.5f, nmi = %.5f, ari = %.5f' % (acc, nmi, ari), ' ; loss=', loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 305,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzsAAAKRCAYAAABz37cTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl4THf///HXZENExK5Cay1tCbVHEUuLUi2qVVXVopZblWpLt1tV/VrV6kJbSt2oWoratShK+VL70sSeiBBLkERkX2Z+f7jNnTTrJJOZyeT5uK65rpM5n3POe+JcrrznfT7vj8FkMpkEAAAAAE7Gxd4BAAAAAEBhINkBAAAA4JRIdgAAAAA4JZIdAAAAAE6JZAcAAACAUyLZAQAAAOCUSHYAIJ0///xTI0eOVNu2bdWwYUPVr19f9evX14IFC+wdWiZvv/22Ob5Lly7ZOxzY0KpVq8z/9qtWrbJ3OADgsNzsHQCAouXatWvavHmz9u7dq+DgYEVFRSkhIUFeXl6qWrWqGjVqpPbt2ysgIEAeHh72Dtcic+bM0fTp0+0dBpzAqlWrFB4eLkkaPXq0naMBgOKLZAdAnty+fVtfffWVli9fruTk5Ez7o6KiFBUVpZMnT2r58uUqX768Ro4cqf79+8vd3d0OEVvm+vXrmjFjhiTJ09NTL7zwgu6//36VKlVKknT//ffbMzwUMatXr9b+/fslkewAgD0ZTCaTyd5BAHBsFy5c0IgRIxQSEmJ+z8/PT23atFH16tXl5eWl6OhohYWFaffu3Tpz5ox53I8//qhWrVrZI2yLrF27VuPHj5ckvfHGGxo2bJidI0JRNnDgQHOyc/r0aTtHAwDFF5UdADmKiorSSy+9pMuXL0uS6tevrw8//FAPP/xwluMnTJig48eP68svv9SePXtsGWqBXL161bz9wAMP2DESAABgLSQ7AHL09ttvmxOdhx9+WD/88IO8vLxyPMbPz0/z58/XggUL5OZWNP6bSf9oXlGbawQAALLGY2wAsnXkyBE999xzkqTSpUtr/fr18vX1tdr5jx07ppUrV2r//v2KiIiQyWRSxYoV1axZM/Xq1Uv+/v45Hl+/fn1JUsuWLbVo0SLFx8dr6dKl2rhxoy5evKiUlBT5+vqqU6dOGjp0qMqWLZvtOXJy9/zSneRv9erVkqRt27apevXq2R6Xl7FJSUlauXKltm7dqrNnzyo6Olru7u4qV66cypUrp/r166tdu3bq3LlzpiTMkliCg4O1dOlS/fXXX7py5YpSUlJUoUIF+fn56YknntBjjz2W4++gU6dOCg8Pl6+vr7Zv367U1FT98ssvWrNmjUJCQpSQkKCqVauqbdu2GjZsmKpWrZrj+XKzatUqvfPOO5KkTz75RH369NHff/+tn376SQcOHNCNGzdUoUIFPfzwwxoxYkSGOVVpaWn69ddftWLFCoWEhOjWrVuqVq2aHnvsMY0YMSLHZD0xMVG7du3Snj17FBgYqLCwMMXGxqpkyZKqWrWqmjdvrv79+6tBgwZZHp/+8bWcvPrqqxnm8vzzXr5165aWLVumLVu26NKlS4qOjlbv3r01derUbH8/d0VGRurJJ5/U9evX5ebmpiVLlqhx48ZZxpGcnKznnntOQUFBkqRp06bpqaeeyjV+ACgqisZXrgDsYuHChebtPn36WC3RSU1N1Ycffqjly5dn2nfx4kVdvHhRa9as0eOPP66pU6eqZMmSuZ7z4sWLGjFihM6dO5fh/XPnzuncuXPauHGjfvzxxxwTAlsLCwvT0KFDdeHChQzvp6SkKD4+XuHh4QoMDDQnFfl9vG7GjBmaPXu20tLSMrx/+fJlXb58WZs2bVLLli01c+ZM+fj45Hq+yMhIjRo1SocPH87w/oULF3ThwgVt2LBB//nPf9SwYcN8xZuVn376SZ988olSU1Mzxf/777/ru+++U7t27RQbG6vXX39df/75Z4bjQ0NDNXfuXG3btk2LFy9W+fLls7xO9+7dzV3U0ouNjTXfS8uWLdPw4cM1btw4q32+9IKCgjRq1ChduXIlX8eXL19eU6dO1dChQ5Wamqo333xTq1evzjLJ+/LLL82JTs+ePUl0ADgdkh0AWTKZTNq7d6/5Z2v+ETR+/Hht3LhRklSiRAn16tVLTZs2lYuLiwIDA7Vy5UrFxcXpt99+0+3bt/XDDz/IYDBke77Y2FgNGzZM58+fV+fOndWuXTuVLVtWly5d0tKlS3X58mWFh4drwoQJWrx4cYZjv/32W0nSxo0b9euvv0qSxowZk6FSkJcEwFImk0ljxowxJzoPPPCAunbtqho1asjNzU0xMTEKDg7Wvn37dPLkyXxfZ/r06ZozZ44kydXVVd27d1fr1q1VsmRJnTlzRr/88otu3Lih/fv368UXX9SKFStUokSJbM+Xmpqq1157TYcPH1arVq306KOPqlKlSrp27ZpWrlyps2fP6tatWxo3bpw2bNhglUcCd+zYoS1btqh8+fJ65plnVK9ePSUmJur333/Xjh07lJycrNdff13btm3ThAkT9Oeff6pp06bq1q2bKlWqpMuXL2vJkiUKDw9XSEiIPv74Y33++edZXispKUk+Pj5q06aNHnjgAVWpUkXu7u66du2agoKCtGnTJqWkpOj7779X+fLl9dJLL2U4fsyYMYqOjtZXX32ls2fPSvrfPZZerVq1srx+dHS0/vWvf+nq1asKCAhQQECAypUrp4iICIt+Z23bttVLL72k+fPnKywsTJMnT9a0adMyjNmzZ4/mz58vSapevbomTZpk0TUAoCgg2QGQpZCQEEVHR0uSSpYsabVJ+7/++qs50alYsaIWLlyounXrmvc/+eSTGjRokF588UVdunRJu3fv1pIlSzRgwIBsz3nixAm5u7tr1qxZ6tixY4Z9zzzzjPr27atLly7p4MGDOn78uPz8/Mz7H330UUnKkFA0a9as0DvIBQYG6sSJE5Kkjh076ttvv5Wrq2uWY8+dO6cKFSpYfI0jR45o7ty5ku60054zZ45atGiRYczgwYM1ZMgQBQYG6vTp0/rqq680YcKEbM957do1Xbt2TZMnT1a/fv0y7Ovfv78GDhyoY8eO6cKFC9q6dau6d+9ucdz/tHnzZvn5+emHH37I8Chi37599e9//1vLly/X7du39fLLLysoKCjLbnq9e/fWU089pevXr+vXX3/V+PHjVbly5UzX+uSTT9SmTZts55q9/vrrGjp0qEJCQjRjxgz17ds3Q8WkefPmkjJWRe/eY3lx5swZubq66quvvtLjjz+e5+OyMm7cOO3bt08nTpzQ2rVr1b59ez3xxBOS7lTnxo8fL5PJJDc3N02fPj3XuXgAUBS52DsAAI7p2rVr5u1q1apZrdHA3T++Jenjjz/OkOjc5evrqy+//NJczZk3b16mR7D+aeTIkZkSHUkqV66cRowYYf55165d+Q3dqsLCwszbTz/9dLaJjiTVrVtX5cqVs/ga8+bN091pmW+99VamREe6U7WaMWOGeT2hZcuWKSYmJsfzPv3005kSHelOlW7s2LHmn3fv3m1xzFlxd3fXV199leWcq1dffdV8nwQFBal9+/ZZtg2vUKGCXnjhBUl35vRk1ymwffv2Od7rvr6++uCDDyRJcXFx2rZtm8WfJzcDBw4scKIj3Wm0MX36dPO/7aRJk3Tp0iVJ0nvvvafr169LkkaNGqUmTZoU+HoA4IhIdgBk6W5VR5K8vb2tcs5Lly6Zqxn333+/AgICsh3r5+en1q1bS5LCw8PN8wqy4urqav5DNit3zyPdmajvCNLPQ7r7uJM1JScna+fOnZLuJDR9+/bNdqyvr6969OghSYqPj881SXnxxRez3de8eXNzsmCt33WnTp2ynS9WpUqVDPtyqgA2a9bMvF2Q2Jo2bWrePn78eL7Pk52c7mVL1a5dW++++66kOwsDv/XWW1q0aJG2b98uSWrRokWGLwMAwNmQ7ACwmfR/GLZt2zbX8Y888oh5+9ixY9mOq1mzZpbf+t9VpUoV8/atW7dyva4tNGvWzJzwfPvtt5o6dapOnTpltfOfOnXK3E67VatWuc6dSf+7zukP+FKlSuXYwc7Dw8NchbLW7zr9Y4dZqVixYp7Gpn8UMKfYbt68qXnz5mnw4MFq3769mjRpovr165tfjRo1Mo9Nvz6TNVSpUkU1atSw6jmfffZZde3aVZJ0+PBhTZkyRZJUtmxZffbZZ3Jx4U8BAM6LOTsAspR+Un5ujzXl1d3HZqQ7CUpu0k/iTn/sP+X2iFf6P/TTr6djTz4+PnrnnXc0adIkpaamav78+Zo/f765nXLz5s3Vvn171alTJ1/nTz+h3Zq/ax8fnxybRUj/+31b63edW4OI9P++OY3Ny33w66+/auLEibp9+3aeYouNjc3TuLxKn5hb00cffaRjx45lSM4mT56se+65p1CuBwCOgmQHQJbST96+fPmyUlNTCzxvJy4uzrx9dx5BTjw9PbM89p+K6jfTzz33nGrXrq3vvvtO+/btk9Fo1M2bN7V161Zt3bpVU6dO1cMPP6x333031+rGPznT79qSaxYkvgMHDuiNN96Q0WiUJD300EPy9/fXvffeqzJlymRIlkaNGiVJ5rHWkpc26/lRunRpVa5c2ZzslC1bVm3atCmUawGAIyHZAZClOnXqyMfHR9HR0UpMTNTJkyczPL6TH6VLlzZvJyQk5Do+Pj4+y2OLirz8IdyyZUu1bNlSUVFROnTokI4cOaIDBw7o77//ltFo1JEjR/T8889r3rx5FnWIK26/a2uYOXOm+d/so48+0rPPPpvluPS/q6Ji5syZGR5PvHXrliZOnKivvvrKjlEBQOErml+HAih0BoNB/v7+5p/Xrl1b4HNWqlTJvB0aGprr+PRjsmoTbA/pv91PSUnJcWxUVFSez1uuXDk9+uijeuutt7R8+XLt2LHD3CY4JSVFn376qUVxpv99FdXftS0lJyfr0KFDkqSGDRtmm+hIdyqdRcmBAwfMay35+vrqoYcekiT99ttvWrVqlT1DA4BCR7IDIFvpu26tWrUqy5XlLZH+UazsWv+m93//939ZHmtPZcqUMW/ntNBjWlqaAgMD832dKlWq6NNPPzUniEFBQUpMTMzz8Q0aNDAnZvv37881MUv/uy5oBa8oio6OVmpqqiTp3nvvzXFsXlpqp5/XdLf9tz3ExMRo/PjxMhqNcnV11bRp0/TFF1+YH1v86KOPzAvbAoAzItkBkK2mTZua20PHxcXpjTfesGhC9oIFC3T48GHzz9WrVzd/q3zq1KkMf2D/099//62//vpLUsZvo+0t/bpAd+PLysaNGxUZGVmga7m5uWWYsH73j/G88PDwUIcOHSTdqTCtXr0627FXrlwxL/Tq6empdu3a5S/gIiz9vKb0ayD9U2xsrBYsWJDr+dLPgbLnY2///ve/zZWo4cOHq3nz5qpZs6bee+89c2xvvPFGrskwABRVJDsAcjR16lRVrVpVkszzR44ePZrjMcePH9fgwYP1ySefZPojaujQoebtCRMmZLneyeXLlzVu3Djz/IkhQ4bkuOimLbVp08Ycy5IlS7Ksdv3999/m9r7ZWbdunX755ZccqzVHjx7VyZMnJUk1atSweIX7IUOGmCfsT5061fyYVnq3bt3Sa6+9Zv6D/LnnnstQvSouypQpY+5aFxgYqN9//z3TmLi4OI0ZM0ZXrlzJ9XzVq1c3b99dW8rWVq5cqU2bNkmSmjRpYm6qIEl9+/ZVt27dJN25X2fOnGmXGAGgsNGgAECOypcvrwULFmjEiBEKDQ3V6dOn1a9fPzVu3Fht2rSRr6+vvLy8dOvWLYWFhWnXrl06c+ZMtufr3r27tm7dqo0bN+r69evq06ePevfurSZNmsjV1VWBgYFauXKluYLUtm1bPf/887b6uLmqUqWKnnjiCa1du1bR0dHq27evnn/+edWpU0fx8fHav3+/Nm7cqLJly6p169bZVn8uXLigb775RlOmTFGbNm3UqFEjVa1aVR4eHoqMjNTBgwe1detWpaWlSbrzrbylmjRpoldeeUXff/+94uLiNHDgQPXo0UOtW7dWyZIldebMGa1cuVI3btyQJNWvX19jxozJ/y+niHvhhRfMSeprr72mnj17qlmzZipdurTOnj2rVatWKSIiQr169dKaNWtyPJe/v78WLVokSXrvvfc0aNAg+fr6mpPP++67T/fdd1+hfZbQ0FD9v//3/yTdaTjx+eefZ+qmeLcd9ZUrVzR37ly1bdtWLVu2LLSYAMAeSHYA5KpWrVpasWKFvvjiC61cuVIpKSk6duxYjgt9VqpUSSNHjsywav1d06ZNk6enp1asWKHExEQtXbpUS5cuzTSua9eumjZtWq7rutjau+++qzNnzujkyZOKjIzUN998k2F/pUqV9O2332b5me66+5ni4+PNraaz4u7urtdee03PPPNMvmIdN26cXF1d9f333ystLU3r1q3TunXrMo1r2bKlZsyYUWitj4uCF154QceOHdP69etlNBq1du3aTI05OnfurA8//DDXZKdDhw5q1qyZDh06pAsXLmjy5MkZ9r/66qsaPXq01T+DdKehxRtvvGGu1k2cODHLhUq9vb312Wef6cUXX5TRaNT48eO1du3aHBfoBYCihmQHQJ54e3tr0qRJGjFihDZt2qS//vpL586dU1RUlBITE+Xl5aVq1aqpUaNGCggIUEBAQLbr8ri5uWnKlCnq27evVqxYoQMHDuj69esyGo2qWLGimjZtqj59+mToBudIfHx8tGzZMv3444/67bffzJ3MqlWrpkcffVSDBg1S+fLlc0x2RowYoZYtW2rv3r06fvy4zp8/r5s3byo1NVWlS5fWfffdp1atWumZZ54pcAVgzJgx6tGjh5YtW6a9e/fqypUrSklJUfny5dW4cWM98cQT6tKlS4Gu4QwMBoM+//xzdejQQcuXL9fJkyeVkJCgChUq6IEHHtCTTz6p7t275+lcrq6umj9/vhYuXKjt27crJCREsbGx5kpdYfr666/NzTF69OihXr16ZTu2RYsWGjZsmGbPnq0rV67o3//+t2bMmFHoMQKArRhM9mwTAwAAAACFhAYFAAAAAJwSyQ4AAAAAp0SyAwAAAMApkewAAAAAcEokOwAAAACcEq2n8ynlRoi9Q0AxUapaO3uHAACA00hNDrd3CNmy5d+X7hVr2+xa9kRlBwAAAIBTorIDAAAAOAJj4S88XNxQ2QEAAADglKjsAAAAAI7AZLR3BE6Hyg4AAAAAp0SyAwAAAMAp8RgbAAAA4AiMPMZmbVR2AAAAADglKjsAAACAAzDRoMDqqOwAAAAAcEpUdgAAAABHwJwdq6OyAwAAAMApUdkBAAAAHAFzdqyOyg4AAAAAp0RlBwAAAHAExjR7R+B0qOwAAAAAcEpUdgAAAABHwJwdq6OyAwAAAMApUdkBAAAAHAHr7FgdlR0AAAAATonKDgAAAOAATMzZsToqOwAAAACcEskOAAAAAKfEY2wAAACAI6BBgdVR2QEAAADglKjsAAAAAI6ABgVWR2UHAAAAgFOisgMAAAA4AmOavSNwOlR2AAAAADglKjsAAACAI2DOjtVR2QEAAADglKjsAAAAAI6AdXasjsoOAAAAAKdEZQcAAABwBMzZsToqOwAAAACcEpUdAAAAwBEwZ8fqqOwAAAAAcEpUdgAAAAAHYDKl2TsEp0NlBwAAAIBTItkBAAAA4JR4jA0AAABwBLSetjoqOwAAAACcEpUdAAAAwBHQetrqqOwAAAAAcEpUdgAAAABHwJwdq6OyAwAAAMApUdkBAAAAHIGRRUWtjcoOAAAAAKdEZQcAAABwBMzZsToqOwAAAACcEpUdAAAAwBGwzo7VUdkBAAAA4JSo7AAAAACOgDk7VkdlBwAAAIBTorIDAAAAOALm7FgdlR0AAAAATolkBwAAAIBT4jE2AAAAwBHwGJvVUdkBAAAA4JSKVGUnISFBR48eVXBwsM6fP6/o6GjFx8crKSlJJUqUkKenp3x8fFSrVi3VqVNHTZo0UalSpewdttNbs/F3vf/xFzmOcXFx0fFdG7PdP/GTr7Rqw2ZJ0q8/z9O91atZNUYUH76+92jSB2+qa5cOqlChnK5cidDadZv10ZQvFB19y97hwclwv8GWuN+cn8mUZu8QnE6RSHZ27dqlJUuWaO/evUpKSsrzcSVKlJC/v78GDBigtm3bFmKExVuDerU1cvCALPcdPhaofYeOqW3r5tkev2P3X1q1YbM8S5VSfEJCYYWJYqB27fu0a+daValSSWvXbdLp0+fUovnDGvPaUHXt2kHtA3opMjLK3mHCSXC/wZa434D8cehk59q1a3rzzTd18OBBSZLJZLLo+MTERO3YsUM7duxQixYt9Nlnn6lKlSqFEWqx1uD+Ompwf50s9w0Y9rok6ZknH89yf2RUtD74dIa6dW6vG5FROnjk70KLE87vmxkfq0qVShoz9n19+9188/ufT/tAY8cO00eTJ2jUq2/bMUI4E+432BL3WzHBnB2rM5gszSBsJDg4WC+//LKuX7+eIckpVaqU6tevr6pVq6pixYoqWbKkPDw8lJycrMTERN28eVNXrlzR6dOnlfCPKkHlypW1YMEC1a5du8DxpdwIKfA5nN2Z4PPq8+K/VKVSBW35ZaFcXV0zjXntnck6FnhKa3+arbHvTdHBI3/zGNs/lKrWzt4hFAm1a9+nM6f26Pz5MN3foE2G/ze8vErrUtgRGQwG3ePrp/h4KogoGO432BL3m3WlJofbO4RsJez4j82uVarDYJtdy54csrKTnJyscePGKSIiQpLk4eGh3r17q3fv3mrUqFGWfzT/U1pamoKCgrRq1SqtXr1aycnJioiI0Lhx47RixQq5u7sX9sco9lau/U2S1PuJrln+m63Z+Lu2/7lXM6ZOlE9Zb1uHByfTIaCNJOn3rX9mqgLHxsZpz54D6tKlg1q3aqbtf+y2R4hwItxvsCXut2LERGXH2hyyG9uGDRt0+vRpGQwG3XvvvVq9erU+/PBDNWnSJE+JjiS5urrKz89PkyZN0qpVq1SjRg1J0unTp7Vhw4bCDB+SEpOStGHLH3J1ddHTPbtl2n/56jVN/Xq2nujaSZ3a+dshQjib+v99lPLs2ayrrmfPnZck1atX8MouwP0GW+J+A/LPIZOdjRvvdO1ycXHRrFmzVKdO1vNB8qpOnTqaNWuWXFzufFySncK3edufirkdq0daNdc9VSpl2Gc0GvXulOnyLFVK74wdYacI4Wy8y5aRJN26FZPl/piY25IkHx+qiCg47jfYEvdbMWI02u5VTDhkshMSEiKDwaAWLVoUONG5q06dOmrVqpVMJpOCg4Otck5kb8W6TZKkZ5/K3Jjgx59X6+CRvzVpwmsq613G1qEBAACgmHDIOTs3b96UJPOjZ9bi6+srSYqMjLTqeZHRuZALOvr3CVWpXFHt/Ftk2Bcadkkz5ixUrx6PqX2blnaKEM4o5tadbzbLZjP/y/u/iXV0dNbfjAKW4H6DLXG/FSPM2bE6h0x2ypYtqxs3bpgbFFjL9evXzedH4Vmx7k5jgj5ZNCYIDg1TcnKK1mz8XWs2/p7l8d37DZEkff3Jv9W5fZvCDRZO4/SZOxXb7J5Zr1e3lqTsn3kHLMH9BlvifgPyzyGTnZo1a+r69evau3evIiIiVLly5QKf89q1a9q7d68MBoPuu+8+K0SJrCQlJWv9pm1ydXVRnye6ZtrvW7VKlu9L0p979+vGzSh17dROpT095VuVNZGQdzt27pEkPfZoexkMhkytWdu0aaG4uHj9te+QvUKEE+F+gy1xvxUjxWguja04ZLLTpUsXHThwQCkpKRo1apTmzJmjcuXK5ft80dHRGj16tJKTk2UwGNStW+buYLCOzX/sUsztWAU80jJTYwLpzgKkk98Zm+WxL706XjduRmnM8JdYZwcWCwm5oC1bdqhLlw7618iXMiy6N2nim/LyKq3v5yxiDQpYBfcbbIn7Dcg/h1xUNDY2Vj179tTVq1clST4+Pho6dKh69eqlChUq5Pk8kZGRWrNmjebNm2eep3PPPfdo/fr1Kl26dIFiZFHRrL048k0dPh6kbz79QB3atrbo2JdeHc+iollgUdG8q137Pu3auVZVqlTS2nWbdOrUWbVs0VQdOz6i02eC1a79U4qMjLJ3mHAS3G+wJe4363HoRUU3f2Oza5Xq+qrNrmVPDpnsSNLBgwf1yiuvKDExUSaTSQaDQZJUq1YtPfjgg6pWrZoqVKigkiVLyt3dXSkpKUpMTNTNmzd1+fJlnTx5UufPn5fJZDKXe0uVKqUffvhBzZo1K3B8JDuZBYeG6akBw1WlckVtWbkgz2si3UWykzWSHctUr15Nkz54U127dFCFCuV05UqE1qzdpI+mfKHo6Fv2Dg9OhvsNtsT9Zh0kO3eQ7DiAoKAgjR49WpcvXza/dzfpyYv0H83X11czZ87Ugw8+aJXYSHZgKyQ7AABYD8nOHcUl2XHIOTt3PfTQQ9q8ebOWLVumpUuXKiQkRJbmZrVr19bzzz+vfv36yd3dvZAiBQAAAAqIBgVW59DJjiS5u7tr4MCBGjhwoEJDQ7Vv3z4FBwcrNDRU0dHRio+PV1JSkkqUKCFPT0+VLVtWtWrVMi8iWrNmTXt/BAAAAAB24PDJTno1a9YkeQEAAIBzorJjdS72DgAAAAAACkORquwAAAAATstEZcfaqOwAAAAAcEpUdgAAAABHwJwdq6OyAwAAAMApUdkBAAAAHAFzdqyOyg4AAAAAp0RlBwAAAHAEzNmxOio7AAAAAJwSlR0AAADAETBnx+qo7AAAAABwSlR2AAAAAEfAnB2ro7IDAAAAwClR2QEAAACQL3v27NGmTZt05MgRRUREKD4+XhUrVlTlypXl5+en1q1by9/fX56enjmeJzY2VqtXr9bmzZsVGhqqW7duqUKFCqpTp466d++unj17ysPDw+L4DCaTyZTfD1ecpdwIsXcIKCZKVWtn7xAAAHAaqcnh9g4hWwnLJ9vsWqWenVig40NCQvTBBx9o//79uY79+uuv1a1bt2z3Hzx4UG+99ZYuX76c7ZggDKIHAAAgAElEQVT69evryy+/VJ06dSyKk8oOAAAAgDw7fvy4Bg8erNu3b0uSDAaDHnzwQdWsWVNeXl6KjY1VSEiIzpw5o7S0tBzPFRQUpFdeeUXx8fGSJHd3d/n7+6ty5cq6ePGiDhw4IKPRqNOnT2vQoEFauXKlqlatmudYSXYAAAAAR1AEHri6cOGChgwZYk50nnjiCb3xxhuqVq1aprExMTHavn27fH19szxXcnKyRo8ebU50HnroIX333XcZkplz585pxIgRunjxoq5fv64333xTP/30U57jpUEBAAAAgFyZTCa9//77iomJkSQNGzZM06dPzzLRkSRvb2/16tVLjRo1ynL/zz//rPDwO48V+vj4aO7cuZmqNnXr1tXs2bPN83UOHDignTt35jlmkh0AAADAERiNtnvlw5YtW8xzdBo1aqTXX3+9QB93yZIl5u0hQ4aoQoUKWY6rW7euevfubf556dKleb4GyQ4AAACAXP3888/m7cGDB8vFJf+pxPnz5xUS8r+GX+mTmayk379nzx7FxcXl6TokOwAAAIAjcODKzvXr17Vnzx5Jd5oIdO7cuUAfdd++febtWrVqqVKlSjmOb9Sokbl9dVJSko4ePZqn69CgAAAAAECOjh07prsr1tSpU0clSpRQXFycVqxYoV9//VUXLlxQYmKiKlSoID8/P/Xo0UOPPfZYtucLDg42bz/00EO5Xt/NzU3333+/OckJDg7WI488kvtxuY4AAAAAUPhM+ZtLYwt///23efuee+7RqVOnNHr0aIWFhWUYFx4ervDwcP32229q2rSpZs6cqYoVK2Y6X2hoaIbz5UX65gXnz5/P0zE8xgYAAAAgR1evXjVv3759W0OHDlVYWJgMBoP8/Pz09NNPq2fPnqpevbp53OHDh9W/f39z97b0oqOjzdtZJUNZSf+o261bt/J0DJUdAAAAwBHks0tafsTExGSZhHh7e8vb2zvL8XcdPHhQ0p3kY8aMGWratKl5n8lk0tKlSzVlyhSlpaUpLCxMkydP1ueff57hfHfX1pGkEiVK5Cnm9OPSH58TKjsAAABAMbNw4UJ17tw502vhwoVZjk9ISMjws5ubm77//vsMiY4kGQwGPf/883rzzTfN723cuFEXLlzIMC4pKcm87e7unqeY7661I0mJiYl5OobKDgAAAOAI/tsAwBYGDRqUZbvnrKo6Uubqy+OPP55jY4EXX3xR8+fPV0REhIxGozZt2qThw4dneb6UlJQ8xZycnGzeLlmyZJ6OIdkBAAAAipnsHlfLzt22z3fl1nrazc1NHTp00PLlyyVJR44cyfZ86as8OUk/7p/xZBtHnkYBAAAAKFw2nLNjKR8fnww/161bN9dj6tSpY96OiIjIsK9s2bLm7Rs3buQphvTj0h+fE+bsAAAAAMhR7dq1M/ycl8pK6dKlzdtxcXEZ9tWqVcu8feXKlTzFkH5c+uNzQmUHAAAAcAQOXNn5ZyUnL93Q0ic4Xl5eGfalr/qcOHEi13OlpqbqzJkzWR6fEyo7AAAAAHLUpEkTlSpVyvzzuXPncj0mODjYvP3PhUNbtWpl3j5//nyuj7IFBQWZE6wSJUqoSZMmeYqbZAcAAABAjkqVKqW2bduaf96+fXuO49PS0rRz507zz82bN8+wv1atWuZH40wmk9asWZPj+VavXm3e9vf3z/CIXE5IdgAAAABHYDLa7pUPAwcONG//9ttvOnnyZLZjFy1apGvXrkm6sz5Ojx49Mo15/vnnzdvz5s1TZGRklucKDg7WqlWrsjwuNyQ7AAAAAHLVqlUrdezYUdKdtXGGDRumo0ePZhhjMpn0888/67PPPjO/N3DgQFWqVCnT+fr16ydfX19JUmRkpIYNG2ZOkO4KDg7WiBEjzG2nW7RooYCAgDzHbDCZbLh6kRNJuRFi7xBQTJSq1s7eIQAA4DRSk8PtHUK24ue8brNreQ77Ml/HRUZGqn///goNDZUkGQwG+fn5qV69ekpOTtaRI0d08eJF8/jmzZtrwYIFcnd3z/J8QUFBGjBggBISEiTdqQL5+/urUqVKCg8P1759+2T8b+OGSpUqacWKFZnm/+SEZCefSHZgKyQ7AABYD8nOHflNdqQ7LaDHjx+v/fv35zjuySef1EcffaSSJUvmOO7gwYN66623dPny5WzH3H///fryyy/ztL5PeiQ7+USyA1sh2QEAwHocOtmZPcZm1/Ic8XWBjjeZTNq2bZvWr1+vwMBA3bhxQ66urqpcubJatmypp59+Wo0bN87z+WJjY7Vq1Spt2rRJYWFhio6OVvny5VW3bl11795dTz75pDw8PCyOk2Qnn0h2YCskOwAAWA/Jzh0FTXaKChYVBQAAABxBPrukIXt0YwMAAADglKjsAAAAAI7AyOwSa6OyAwAAAMApUdkBAAAAHIGROTvWRmUHAAAAgFOisgMAAAA4Aio7VkdlBwAAAIBTorIDAAAAOAIT3disjcoOAAAAAKdEsgMAAADAKfEYGwAAAOAIaFBgdVR2AAAAADglKjsAAACAIzDSoMDaqOwAAAAAcEpUdgAAAABHYGLOjrVR2QEAAADglKjsAAAAAI6AOTtWR7KTT2lhgfYOAQAAAEAOSHYAAAAAB2BinR2rY84OAAAAAKdEZQcAAABwBMzZsToqOwAAAACcEpUdAAAAwBGwzo7VUdkBAAAA4JSo7AAAAACOgDk7VkdlBwAAAIBTItkBAAAA4JR4jA0AAABwBCwqanVUdgAAAAA4JSo7AAAAgCOgQYHVUdkBAAAA4JSo7AAAAACOgEVFrY7KDgAAAACnRGUHAAAAcATM2bE6KjsAAAAAnBKVHQAAAMABmFhnx+qo7AAAAABwSlR2AAAAAEfAnB2ro7IDAAAAwClR2QEAAAAcAZUdq6OyAwAAAMApUdkBAAAAHIGJbmzWRmUHAAAAgFMi2QEAAADglHiMDQAAAHAENCiwOio7AAAAAJwSlR0AAADAAZio7FgdlR0AAAAATonKDgAAAOAIqOxYHZUdAAAAAE6Jyg4AAADgCIwsKmptVHYAAAAAOCUqOwAAAIAjYM6O1VHZAQAAAOCUqOwAAAAAjoDKjtVR2QEAAADglKjsAAAAAA7AZKKyY21UdgAAAAA4JSo7AAAAgCNgzo7VUdkBAAAA4JRIdgAAAAA4JR5jAwAAABwBj7FZHZUdAAAAAE6Jyg4AAADgAExUdqyOyg4AAAAAp0RlBwAAAHAEVHasjsoOAAAAAKdEZQcAAABwBEZ7B+B8qOwAAAAAcEpUdgAAAAAHQDc266OyAwAAAMApUdkBAAAAHAGVHaujsgMAAADAKVHZAQAAABwB3disjsoOAAAAAKdEZQcAAABwAHRjsz4qOwAAAACcEskOAAAAAKfEY2ywmn2BZ7V08//p+NkLiolLkI9XadW9t6oGdGurdg8/YB6XnJKqVX/s07o/Dyn82k0lpaSqagUftW5UTy/2CFC1SuXs+ClQ1Pn63qNJH7yprl06qEKFcrpyJUJr123WR1O+UHT0LXuHByfD/QZb4n4rBmhQYHUGk8nEw4H5kHh4nb1DcChfLt6gBRt2qkr5smrbpIF8ypRWVEysTpwPV+uGdfX6gCckSalpaRry0WwdPR2qWtUqq1XDevJwd1VQ8CUdOhWiMp4ltfDDV1WnehU7fyLH4dV6pL1DKDJq175Pu3auVZUqlbR23SadPn1OLZo/rI4dH9Gp0+fUPqCXIiOj7B0mnAT3G2yJ+816UpPD7R1CtqKe7mCza5X7ZYfNrmVPxbKyM3z4cJ07d04Gg0Fbt261dzhF3i/b9mnBhp16sn0zTXylr9zdMt5WKalp5u3tBwJ19HSoWjWsq9nvvCIXl/89Sfndis36ftVWLdywU5NHPGuz+OE8vpnxsapUqaQxY9/Xt9/NN7//+bQPNHbsMH00eYJGvfq2HSOEM+F+gy1xvxUPNCiwvmI5Z+fatWsKDw9XeLjjZvZFRXJKqmYu/033VPTJMtGRJHc3V/P2pYhISVK7hx/IkOhIUsfmD0mSom7HFmLEcFa1a9+nLl066Pz5MH03a0GGfZMmf67Y2Di9MOBpeXqWsk+AcCrcb7Al7jcg/4plsgPr2fv3GUXFxKlzi0YyGAz68/BJ/WfdH1r82y4dOxOaafzdx9N2Hz0lozHjg6k7D5+UJLVuWK/Q44bz6RDQRpL0+9Y/9c+nc2Nj47RnzwGVLu2p1q2a2SM8OBnuN9gS91sxYrThq5golo+xwXqCgi9Kkjzc3dTvna907uLVDPubNaitz18fqPLeXpKk9g8/oM4tG2rb/kA9Pf4LtW5UT+5urjoRcklHToeqf9dH1K9LG5t/DhR99e+vI0k6ezYky/1nz51Xly4dVK9ebW3/Y7ctQ4MT4n6DLXG/oagYPXq0tmzZYv65ZcuWWrRoUbbjZ86cqW+++caiayxevFjNmzfP83iSHRRIZMydR84Wbtip2r6VNf+Df6lBzWoKj4jU9MUbtPf4Gb311SLNm3hnkr3BYND0sS9q9i+/a+7qbQoJv2Y+V6uGddX9kYfl5uqa5bWAnHiXLSNJunUrJsv9MTG3JUk+Pt42iwnOi/sNtsT9VnyYinDFZevWrRkSHUfh0MmOpZleXt24caNQzlscGf9bTnd1ddHXb70s30rlJUn17r1HX44bpKfGTdPBkyE6diZUje+vqaTkFL3/3TLtPnZa77zcSx2bP6SSHh46eiZUny5cq5c/nKXPx76gjs0b2vNjAQAAII9iY2M1efLkAp2jUaNG8vPzy3VclSqWdex1+GTHYDDYOwzkoMx/J0M2qFnNnOjcVaqEh9o0rq/Vf+xXYPBFNb6/pv6z7g9t2Xdc4wc9pWce9TePbdukgSqV89azb3+paQvXkezAYjG37nyzWbZs1t9senvf+WY0Ojrrb0YBS3C/wZa434qRIlrZmT59uq5duyZ3d3d16tRJmzdvtvgcAQEBGj16tNVjKxINCkwmk1VfsJ6a91SS9L+k55+8S995PzE5RZL053+bELR8sE6msfXvqybv0qV0+UaUom/HFUa4cGKnzwRLkurVq53l/np1a0nK/pl3wBLcb7Al7jc4ssOHD2vp0qWSpJdffln16jlWoymHruzcreoYDAY1btxY7u7uVjlvYGCgEhISrHKu4q5Vw3oyGAwKCb8mo9GYqZ303YYFvpXvVH2SU1MlSZExmZOZ5JRUxScmScrYrhrIix0790iSHnu0vQwGQ4YvNry8SqtNmxaKi4vXX/sO2StEOBHuN9gS91vxUdTm7KSkpGjixIkymUyqXr26Ro0apblz59o7rAwcOtmpWbOmzp8/L4PBoLfeekvNmlmnpWKvXr106tQpq5yruKtWqZwCmj6gHYdOaPGm3RrYvb15357jp7Xn+BmVKV1KjzSuL0lqWr+Wzl28qnlrt+nh+jXl4f6/W3DWyi1KTTPqoTo1VLpUSZt/FhRtISEXtGXLDnXp0kH/GvlShkX3Jk18U15epfX9nEWKj+eLDhQc9xtsifsNjmrOnDk6e/asJGnixIkqWdLx/n5z6GTHz89P58+flyQdP37caskOrOvdl3vrVOhlfb5ovXYdOakGNX0VHhGpPw4GydXFoEmv9DU/5ja0d2ftPHxC+wLP6ak3pumRxvVVwsNdR0+HKjD4okp6uGvCi0/a+ROhqHr1tXe1a+daff3VFHXq1FanTp1VyxZN1bHjIzp9Jlj/nvipvUOEE+F+gy1xvxUTRaiyExISotmzZ0uSunbtqoCAADtHlDWHT3bWrl0rSTp27Jido0F2qlTw0dKPx+j7X7Zq5+EgHTp5Xl6lSqh90wc05KlOalT33v+NLV9Wyz4Zq/nr/tCuI6e0dudBGY0mVSpXRk8GNNfgnh1Vy7eyHT8NirKQkAtq5d9dkz54U127dNDj3TrpypUIfT3jB3005QtFR9+yd4hwItxvsCXuNzgSk8mkiRMnKjk5WV5eXnrvvfcKfM4bN25ozZo1Cg0NVUJCgry9vVW9enW1aNFC1apVy/d5HT7Zuev48eN2jAS5Ke/tpXde7qV3Xu6Vp7FvvNBTb7zQ0waRobi5dOmyhr4yzt5hoJjgfoMtcb85v6IyZ2fFihU6cOCAJGnMmDEWt4POyrJly7Rs2bIs9/n7+2vMmDF6+OGHLT6vQ3dja9CggTw8PGQymXTlyhXdvHnTKuelIxsAAABguevXr+uzzz6TJD300EMaMGBAoV9z7969GjBggBYsWGDxsdlWdqKjowsSUwY+Pj75Os7d3V2TJ0/WrVt3SrOp/+3kVVB3H40DAAAAiqOYmBjFxGRem8nb21ve3lmv6SRJU6ZMUUxMjFxcXDR58mS5uhasg27dunXVrVs3tWzZUnXr1pW3t7cSExN1/vx5bd26VYsXL1ZsbKzS0tL0ySefqFy5cnrqqafyfP5skx1/f//sdlnEYDDoxIkT+T6+V6/cH4sCAAAAijpbPsa2cOFCffPNN5nef/XVV7Nd3HP79u3atGmTJGnAgAFq2LBgi8APHDgwy2u5u7vLz89Pfn5+euaZZzR06FCFhoZKupNsdezYMceELL1skx0e9QIAAACc06BBg9S7d+9M72eXRMTFxWny5MmSpMqVK2vs2LEFjiEvT3/VqFFDs2fPVs+ePZWSkqKYmBgtX75cQ4cOzdM1sk121q9fn/dIAQAAABSILSs7uT2u9k9ffvmlrly5Ikl677335OXlVVihZVKrVi316NFDa9askSTt2rWr4MlOvXr1rBMdAAAAgCIrKChIixcvliQFBASoW7duNo/B39/fnOwEBwfn+TiHbj0NAAAAFBsmg70jyNLp06dlNN4pO124cEHPPvtstmOvXr1q3g4KCsowdtq0aapZs2a+YqhYsaJ525JGagVOdpKSkpSQkJDvjmsAAAAAioa7jQLyIi4uTseOHTP/nJCQkO/rJiYmmrc9PT3zfFy+kp3g4GB9//332r17t6KiojJ1XLt165ZmzZolg8Gg119/XR4eHvm5DAAAAFBsFJVFRe0hfa5RqVKlPB9ncbKzadMmTZgwQcnJydl2bCtbtqwCAwN16NAhNW7c2C7P9QEAAAAouD59+qhPnz55Gjtz5kxzS+uWLVtq0aJFBb5+ampqhuZpzZs3z/OxLpZcKDQ0VBMmTFBSUpL69OmjRYsWZfv4Wu/evWUymbRz505LLgEAAAAUSyajwWYve4uLi8vz2OnTpyssLMz8c8+ePfN8rEWVnf/85z9KSkpS3759NWXKlDsncMv6FG3atJEkBQYGWnIJAAAAAE5u/vz5OnLkiPr376/27dtnOe3l6tWrmj59utatW2d+79FHH7WosmNRsvPXX3/JYDBo5MiRuY695557VKJECXM/bgAAAADZK05zdkwmk3bv3q3du3erZMmSql+/vmrUqKHSpUsrKSlJoaGhCgwMVGpqqvmYBg0a6NNPP7XoOhYlO9euXVPJkiXl6+ubp/GlSpVSbGysRQEBAAAAKD4SExN17NixDJ3b0nN1ddXTTz+td955x6JObJKFyY6bm5uSk5PzNDY5OVmxsbE2XV0VAAAAKKpMDrrOTmEYMmSImjdvrsOHD+vYsWO6cuWKoqKiFB0dLTc3N3l7e6tu3bpq2rSpevfunediyz9ZlOz4+vrq7NmzCgsL07333pvj2L179yo1NVV16tTJV2AAAAAAipbRo0dr9OjRuY7z9PSUv7+//P39CzUei7qxtWvXTpJybSGXkJCg6dOny2AwqEOHDvkODgAAACguTEbbvYoLi5Kdl156SaVKldLixYs1Z84cJSUlZdhvNBq1e/du9e/fX2fOnFG5cuXUv39/qwYMAAAAAHlhMGW3Mmg2/vjjD40ePVppaWkqUaKEUlNTlZaWpnvvvVcRERFKTEyUyWRSyZIl9f3336tVq1aFFbtdJR5el/sgwAq8Wufe/RAAAORNanK4vUPI1qVWnWx2rer7ttvsWvZkUWVHkjp27KglS5aoUaNGSkxMVGpqqkwmky5cuKCEhASZTCY1atRIP/30k9MmOgAAAIC1FadFRW3FogYFd/n5+Wn58uU6e/asDh8+rIiICKWlpalSpUpq2rSpHnjgAWvHCQAAAAAWyVeyc1e9evVUr149a8UCAAAAFFuWTS5BXlj8GBsAAAAAFAX5ruzExMRo27ZtCgwM1M2bNyVJFSpUUMOGDdW5c2d5e3tbLUgAAADA2RWnuTS2YnGyk5aWpu+++07z5s3L1Hr6rhIlSmjw4MEaNWqUXF1dCxwkAAAAAFjK4mRn7Nix2rp1q0wmk9zc3FSvXj1VqVJFknTt2jWdPXtWiYmJmjVrlk6fPq1vv/3W6kEDAAAAzobKjvVZlOysWLFCv//+uyRp0KBBGj58uMqXL59hTFRUlObMmaP58+dr+/btWrFihZ555hnrRQwAAAAAeWBRg4Lly5fLYDBo5MiReueddzIlOpJUrlw5TZgwQSNHjpTJZNLy5cutFiwAAADgrEwm272KC4uSnXPnzsnFxUVDhgzJdeyQIUPk4uKic+fO5Ts4AAAAAMgvix5jc3d3l4eHh7y8vHId6+XlpTJlyshUnFJHAAAAIJ+Ys2N9FlV2ateurdu3bysmJibXsTExMbp9+7bq1q2b7+AAAAAAIL8sSnaee+45GY1GzZ07N9exc+fOldFoVL9+/fIdHAAAAFBcmEwGm72KC4seY+vVq5eOHz+uH374QfHx8Ro+fLgqV66cYcyNGzc0e/ZsLV68WAMGDNBTTz1l1YABAAAAIC+yTXZGjBiR7UFeXl5asmSJli1bplq1amVYZ+f8+fMyGo0qU6aMwsPDNXLkSM2aNcv6kQMAAABOxGS0dwTOJ9tkZ8eOHbkenJaWpnPnzmXZcS0mJkY7duyQwVB8ymQAAAAAHEe2yc7LL79syzgAAAAAwKqyTXYmTJhgyzgAAACAYs1YjBoH2IpF3dgAAAAAoKiwqBsbAAAAgMJRnFpC2wqVHQAAAABOKV+VnfPnz2vZsmU6dOiQrl69qoSEBJlMpizHGgwGHTp0qEBBAgAAAM7OZKSyY20WJzsrV67Uhx9+qNTU1GwTnPRoPQ0AAADAHixKdoKCgvTBBx8oLS1NvXv3VkBAgMaOHauyZcvq448/1o0bN7Rnzx5t3bpV3t7emjBhgnx8fAordgAAAMBp5KGOAAtZlOwsXLhQaWlpeu655zRp0iTz++7u7urcubMkqV+/fjp69KiGDx+uuXPn6pdffrFqwAAAAACQFxY1KDh48KAMBoMGDx6c47gmTZro3XffVXBwsObNm1egAAEAAIDiwGQ02OxVXFiU7Ny4cUPu7u669957/3cCFxclJydnGvv444/Lzc1NmzdvLniUAAAAAGAhix5jK1mypFxcMuZHnp6eun37tpKTk+Xh4WF+38PDQ56enrp06ZJ1IgUAAACcmJF1dqzOospO5cqVFRsbq7S0NPN7NWrUkCSdOHEiw9gbN24oJiYmTx3bAAAAAMDaLEp2ateurbS0NJ09e9b8XvPmzWUymTR79mylpqZKktLS0vTpp59KkurUqWPFcAEAAADnZDIZbPYqLixKdtq0aSOTyaSdO3ea33v++efl5uamnTt3qlOnThoyZIg6d+6sDRs2yGAw6LnnnrN60AAAAACQG4vm7HTr1k0XLlyQp6en+b1atWrp448/1sSJExUREaGIiAjzvhdeeEF9+/a1XrQAAACAk2L2h/UZTFaaVBMREaFt27bp6tWr8vLyUrt27dSgQQNrnNohJR5eZ+8QUEx4tR5p7xAAAHAaqcnh9g4hW8dr9rTZtfxC19vsWvZkUWUnJ5UrV1b//v2tdToAAAAAKBCrJTsAAAAA8o/W09ZnUYMCAAAAACgqsq3sTJkyxWoXef/99612LgAAAMAZFaeW0LaSbbLz008/yWAo2C/cZDLJYDCQ7AAAAACwuWyTnS5dutgyDgAAAKBYo/W09WWb7MyYMcOWcQAAAACAVdGNDQAAAHAAdGOzPpKdfDJFXbV3CAAAAAByQLIDAAAAOAC6sVkf6+wAAAAAcEpUdgAAAAAHwJwd66OyAwAAAMApUdkBAAAAHADL7FgflR0AAAAATonKDgAAAOAAmLNjfVR2AAAAADilfFd2wsLCtHv3bl2+fFmJiYl6//33zfuMRqPi4+MlSV5eXgWPEgAAAHByrLNjfRYnOwkJCZo0aZLWr18vk+l/06jSJzvx8fHq1KmTYmNjtWHDBtWuXds60QIAAABAHln0GFtaWpqGDx+udevWydXVVS1atJC7u3umcV5eXurdu7eMRqM2bdpktWABAAAAIK8sSnZWr16t/fv365577tGvv/6qH3/8Ud7e3lmO7dq1qyTpwIEDBY8SAAAAcHJGG76KC4uSnXXr1slgMOjtt99WjRo1chzbsGFDubi4KDg4uEABAgAAAEB+WDRn5/Tp0zIYDOrQoUOuYz08PFSmTBlFRUXlNzYAAACg2DCJBgXWZlFlJy4uTqVLl5aHh0eexqempsrNjaV8AAAAANieRclOuXLlFBsbq4SEhFzHXr58WXFxcapYsWK+gwMAAACKC6PJdq/iwqJkp1GjRpKkXbt25Tp22bJlkqRmzZrlIywAAAAAKBiLkp3evXvLZDLp66+/VmRkZLbj1q1bp3nz5slgMOjpp58ucJAAAACAszPKYLNXcWHRhJrHHntMHTp00I4dO9S3b1/17t1bSUlJkqT169fr8uXL+uOPP3Ts2DGZTCb17NlTLVq0KJTAAQAAACAnBpPJZNFTe/Hx8Xrrrbe0bds2GQyZs8K7p+vWrZumTZuW52YGRU3Ctjn2DgHFRJnHP7R3CAAAOI3U5HB7h5CtbVX62exana/9bLNr2ZPFrdI8PT317bffaufOnfrll190+PBh3bx5UyaTST4+PmratKn69eungICAwogXAAAAAPIk332hA30dEnAAACAASURBVAICzAmN0WiU0WikzTQAAACQT0Z7B+CErJKduLi4yMXFol4HAAAAAFCoKMUAAAAADsBUjLqk2YpFyc6WLVvydZEuXbrk6zgAAAAAyC+Lkp3XXnstyw5sOTEYDDpx4oRFxwAAAADFDXN2rM+iZKds2bI5JjtxcXFKSUmRJHl4eMjT07Ng0QEAAABAPlmU7Ozbty/XMSdOnNCsWbO0e/duffzxx+rUqVO+gwMAAACA/LJ6C7UHH3xQM2fOVOfOnTVu3DidO3fO2pcAAAAAnI7Rhq/iotD6RY8ZM0aJiYmaNWtWYV0CAAAAALJVaK2na9SooTJlyuTp0TcAAACguKP1tPUVWrKTnJys+Ph4JSYmFtYlAAAAACBbhZbs/P/27js+qir///h7SCOVQBIgFCEFULOIdIIo6lIEFpSVBRekiQKuRNyfoOIKwiL4tWJBURSRVbHAUl0pooJgMNLBUEIagSQIIaQHAsz8/shmNjFtApOZYfJ6+sjjcefec+79DLmPOJ/5nHvO+vXrdeXKFTVt2rS2LgEAAAA4DSOFHaurUbKTlZVV5fGLFy/q9OnT2rRpk5YvXy6DwcBsbAAAAADsokbJTmRkpMVtTSaTWrVqpb/97W81DgoAAACoa4w8s2N1NUp2TCaTRe0CAwM1aNAgPfbYY/Lz87uqwAAAAADgWtQo2Vm/fn2Vx11cXNSgQQMFBARcU1AAAABAXWNZWQE1UaNkp02bNrUVBwAAAABYVY2SnRdeeEGSNG7cOLVo0aJWAgIAAADqIqO9A3BCNUp2li9fLhcXFz377LO1FQ8AAAAAB5SZmam9e/fq4MGDiouLU0pKis6cOaOCggK5u7urQYMGateunXr27KkhQ4aoUaNGFp87Ly9Pq1ev1qZNm5ScnKzs7GwFBAQoLCxMAwcO1ODBg+Xu7l7jmA0mS2cdkNSrVy9dunRJMTExNb6Qsyn8brG9Q0Ad4Ttgjr1DAADAaVwuSrV3CJVaGTzKZtcalv5ZjftMmjRJW7dutaitl5eXnnzyST344IPVtt29e7emT5+utLS0Stu0a9dOCxYsUFhYmKXhSqphZSciIkI//vijzp49q6CgoBpdCAAAAIBzKKm6BAcHy9PTU4WFhUpMTFRsbKyMRqMKCgo0d+5cZWRk6Iknnqj0PLGxsXrkkUdUUFAgSXJzc1NkZKQaN26skydPateuXTIajTp27JjGjh2rlStXqmnTphbHWaNkZ9SoUdq2bZsWLVqkWbNm1aQrAAAAgCo4+mxs3bp10913362ePXuqZcuWFbY5deqUnn/+ee3YsUOS9N577+muu+5Shw4dyrUtKipSVFSUOdGJiIjQu+++WyaZiY+P1+TJk3Xy5EmdPXtW06ZN06effmpxzPVq8gbvuOMO/f3vf9cXX3yhGTNmKCkpqSbdAQAAAFynJkyYoBEjRlSa6EhSixYttGjRIvNwM5PJpBUrVlTY9ssvv1RqavGwQn9/f33wwQflqjbh4eF67733zM/r7Nq1S9u2bbM45hpVdgYPHixJ8vDw0Jo1a7RmzRrzujr16lWcNxkMBq1bt64mlwEAAABwnXJ3d9fgwYP1xhtvSJIOHz5cYbvly5ebtydMmFDpWp3h4eEaOnSovvzyS0nS559/rt69e1sUS42SnePHj5fbl5WVpaysrEr7GAyGmlwCAAAAqJOcaerphg0bmrfz8/PLHU9KSlJiYqL59dChQ6s8X+lkJzo6Wvn5+fL29q42jholOzNnzqxJcwAAAAB1UEJCgnm7ovU5S8/uHBISUu3kZ+3bt5eXl5cKCgp08eJF7d+/X7fddlu1cdR4ggIAAAAA1md0kgFR8fHxWrlypfl13759y7UpnQxFRERUe05XV1e1bdtW+/fvN/e/5mRnzZo18vDw0IABA6o9EQAAAIC6qaCgQCkpKdqyZYuWLl1qnmGtU6dOGjZsWLn2ycnJ5u3g4GCLrlF68gJLJ0qrMtl55plnFBQURLIDAAAA1DKjrp/STnR0tMaPH1/pcYPBoEGDBmnevHlydS2fcpR+5j8wMNCia5Ye6padnW1Rn2qHsZlMjj7jNwAAAICayMnJUU5OTrn9fn5+8vPzu6ZzBwYGav78+VXOmFZS+ZGKZ3q2ROl2pftXpUbP7AAAAACoHbYsMSxbtkwLFy4st3/KlCmKioqqtn/Tpk3Nz/MbjUbl5uYqLi5OcXFxysjI0KRJk9SnTx/Nnj27wsrNxYsXzdtubm4WxVyy1o4kXbhwwaI+JDsAAABAHTN27NgKp3u2tKoTGhqqWbNmldt//PhxzZ49W7t379a3336rY8eO6Ysvvii3hk7pKs2lS5csumZRUZF5u379+hb1qXglUAAAAAA2ZTTY7sfPz08tWrQo93OtQ9jatGmjpUuXqmPHjpKklJQUzZ8/v1w7Ly8v83bpKk9VSrcr3b8qJDsAAAAArMbd3V1PPfWU+fWGDRuUmZlZpk2DBg3M2xkZGRadt3S70v2rUu0wtnPnzummm26y6GQVMRgMOnz48FX3BwAAAOoCo70DsKKOHTvK09NThYWFunLlig4dOlRmwoKQkBBt375dkpSenm7ROUu3CwkJsaiPRc/sMCMbAAAAAEsZDAb5+PiosLBQksrN/BYWFmbetqQwcvnyZcXFxVXYvyrVJjuenp566KGHLDoZAAAAgKvjTOWFS5culVlL5/fDzrp3727eTkpKUkZGRpXr7cTGxpqnm/bw8NCtt95qURzVJjteXl6aMmWKRScDAAAAgOjo6DKzrIWGhpY5HhISotDQUCUmJspkMmnNmjV6+OGHKz3f6tWrzduRkZHy9va2KA6mnobVxBw9oS+27dfBpDTlFFyUv3d9hTcL1Mi7Oun2PxTf4KnnsjVo5oeVnqN/53Z6acKfbBUynFDz5sGa/fw09e93pwICGio9/YzWrtukuS+8rqwsy1ZbBizF/QZb4n5zfkaDvSOoXFZWlvz9/S1qm5ubq1deecX8+uabb1aLFi3KtRs5cqReeOEFSdKSJUv05z//WY0aNSrXLiEhQatWrSrTz1IkO7CKBau2admW3Wri76Pe7cPk7+Op83mFOpLym3bHnTQnOyXatgjSXbeElztPeLPKy5dAdUJDW2n7trVq0iRIa9dt1LFj8erapaOmPv6w+ve/U3f0vk+ZmeftHSacBPcbbIn7Dfa2Zs0affPNNxo9erT++Mc/Vjj1s8lk0s6dOzVv3jzFx8eb9z/55JMVnnPEiBFaunSpUlNTlZmZqYkTJ+qdd95RkyZNzG0SEhI0efJk87TTXbt2LTPRQXWum2QnLy9PP//8s06dOiWj0aimTZuqR48eFWZ/VVm3bp127twpg8FQ4ZzfqLl/7zioZVt2a3CPCM0a2Vduri5ljl+6cqVcn3YtGuvRP/W0VYioIxa+NV9NmgRp6hPP6Z13l5r3v/ry83riiYma+8+n9diUZ+wYIZwJ9xtsifsNjuDAgQM6cOCA3NzcFB4ertatW8vPz09XrlzRuXPn9Ouvv+rs2bNl+jz55JPq1atXhedzd3fX22+/rVGjRqmwsFCHDh1Snz59FBkZqaCgIKWmpiomJkZGY/E8dUFBQWUqRpYwmKqYau3GG29UYGCgduzYUaOTWlNRUZHeeOMNffrpp+VWVzUYDLrnnns0depUtWrVyqLzzZ07V5999pkMBoOOHDly1XEVfrf4qvs6k6JLl9X/H4tV391V62ZPKJfo/F7JMLbBPSI0d8w9Nory+uY7YI69Q7guhIa2UtzRaCUlpajtjT3LzCLp4+OtUyn7ZDAYFNz8FhUUFNoxUjgD7jfYEvebdV0uSrV3CJX6oMWDNrvWI6c+rVH75cuXa84cyz+TNG3aVDNnzlSfPn2qbbt7925Nnz5daWlplbZp27atFixYoPDw8iODquLQlZ3CwkI98sgj2rNnT4XTX5tMJm3YsEE//PCDnn76aT3wwAN2iLJu+/noCZ3PK9SouzvJYDDox0OJSkjPkLurq/7Quqk6hDarsN/ZrDyt3H5AWfkX5O9dX7eENFPbFkE2jh7O5M7exZXCb7f8WO7vRV5evqKjd6lfvzvVo3tnff+D/b7AgXPgfoMtcb/BEYwcOVI9evRQdHS09u/fr4SEBKWlpSk/P1+S5Ovrq+DgYEVERKh3796688475epqWarRpUsXrV+/XqtWrdLGjRuVkpKirKwsNWrUSOHh4Ro4cKCGDBkid3f3GsddZQRHjx6t8Qmtaf78+dq9e7cMhuKntdzc3NSmTRtJ0okTJ8z/uIWFhZozZ4727dunefPmWfwPi2v364nTkiQPV1c98OInik8ruwJu5/AWeuWRwWrkW3Zc589HT+jnoyfK7OvSpqXmjr1HwY38ajdoOKV2bYvn2z9+PLHC48fjk9Sv351q0yaUDwO4ZtxvsCXut7rD0RcVDQ0NVWhoqB580PoVKB8fH40ZM0Zjxoyx6nkdNitISEjQihUrZDAYZDKZNHbsWEVFRcnHx0dS8fC2TZs2acGCBUpLS5PJZNK6deuUmZmpt99+W/Xr17fzO6gbzucWl8uXbdml0KYBWvr/Rqhdi8ZKPZet11dt084jJzT9w/Va8vcRkiRPdzdNHNBDd3UIV/PA4vnWj6dm6L3/RGtX3ElNfHOFvnp2jDw93Oz2nnB98mvgK0nKzs6p8HhOTq4kyd+fZBrXjvsNtsT9Bly9evYOoDL//ve/zdsPP/ywZsyYYU50pOIHmgYPHqwNGzboL3/5i3n/jh07NH78eOXm5to03rrK+N9yuku9enrj0fvUMbyFvOq7q03zIL0+6V418ffRnuOndCCxeAxmI18v/W3wbbrphiby86ovP6/66tymhRZFDVP71sE6eTZLq346ZM+3BAAAYBcmg+1+6gqHTXZ2794tSfLz86tyUVMPDw/NnTtXc+fOlYtL8cPx+/fv1+jRo5WZmWmTWOsyX08PSdKNLRureUDZlXE93d3U8+bWkqRfk09XeR5Xl3oaetsfJEl7409ZP1A4vZzs4i84GjSo+JtNP7/ib0azsir+ZhSoCe432BL3G3D1HDbZOXnypAwGg7p06WLRkLS//OUveu+998xtjx49qlGjRum3336r7VDrtNZNGkr6X9Lze75exb+Pi7+bSa8iDX2Kn+spLKq+LfB7x+ISJElt2oRWeLxNeIikyse8AzXB/QZb4n6rO4w2/KkrHDbZKRmGFhho+SKTvXr10scffyw/Pz8ZDAYlJSXpr3/9q1JSUmorzDqvW7tWMhikxNPnZDSWnzEv4b8TFvy+6lORQ0npFrcFfm/rtmhJUt8+d5gnNSnh4+Otnj27Kj+/QD/H7LFHeHAy3G+wJe434Oo5bLLj4VFcKajpszcdOnTQJ598ooCAABkMBqWlpWnkyJE6fvx4bYRZ5zUL8FPv9mFKz8zVZz/sLXMs+nCyoo8ky9fTQz0jir91OpLyW4VJUczRE/r0++I/0oO63VT7gcPpJCae0ObNWxUScoP+9ui4Msdmz5omHx9vffrZv1mDAlbB/QZb4n6rO6jsWF+Vi4ra04ABA5SUlKSbb75Zq1atqnH/EydOaPz48UpPT5fJZJK/v78WL16stWvXsqiolf12PldjX/1cp8/nqnu7G9SuZWOlncvWDwfiZZBB/zdhkPp0bCtJmrDgS6WcyVKH0GZq0rB4wonjqRn65Vhx9e2xwbfpkQE97PZeHBGLilouNLSVtm9bqyZNgrR23UYdPXpc3bp20l133aZjcQm6/Y57lZl53t5hwklwv8GWuN+sx5EXFV3Y0naLik45WbNFRa9XDpvsPP7449q8ebNcXFz0888/y9fXt8bnSE9P17hx43TiRPF6Lt7e3goPD9eBAwdIdqwsM7dAi7/ZqW2HEnQ2O18+9d3VMbyFHurfTe1bB5vbrf7pkL4/EK/4tAxl5Rfq8hWjAny9dEtIsB64s6M6hbew47twTCQ7NdOiRTPNfn6a+ve7UwEBDZWefkZr1m7U3BdeV1ZWtr3Dg5PhfoMtcb9ZhyMnO2/bMNmJItmxr2XLlunFF1+UwWDQnDlzNHz48Ks6z7lz5zR+/HjFxcVJknndHpIdXC9IdgAAsB6SnWJ1Jdlx2Gd2evXqZd7+8ssvr/o8AQEB+vTTT9WhQwdrhAUAAADUCqPBdj91hcMmO2FhYQoLC5PJZNLhw4f1/fffX/W5/Pz89PHHH6tHjx5y0EIWAAAAACtztXcAVfnHP/6hX3/9VZJkNF7bvBGenp5avHix3nzzTRYbBQAAgMOpS7Ok2YpDJzs9e/ZUz549rXY+d3d3TZ8+3WrnAwAAAOC4HHYYGwAAAABcC4eu7AAAAAB1BcPYrI/KDgAAAACnRGUHAAAAcADMGWx9VHYAAAAAOCUqOwAAAIADqEuLfdoKlR0AAAAATonKDgAAAOAAmI3N+qjsAAAAAHBKVHYAAAAAB8BsbNZHZQcAAACAU6KyAwAAADgAI7Udq6OyAwAAAMApUdkBAAAAHACzsVkflR0AAAAATonKDgAAAOAAeGLH+qjsAAAAAHBKJDsAAAAAnBLD2AAAAAAHwAQF1kdlBwAAAIBTorIDAAAAOACjwd4ROB8qOwAAAACcEpUdAAAAwAEYmXza6qjsAAAAAHBKVHYAAAAAB0Bdx/qo7AAAAABwSlR2AAAAAAfAOjvWR2UHAAAAgFOisgMAAAA4AGZjsz4qOwAAAACcEpUdAAAAwAFQ17E+KjsAAAAAnBKVHQAAAMABMBub9VHZAQAAAOCUSHYAAAAAOCWGsQEAAAAOgKmnrY/KDgAAAACnRGUHAAAAcADUdayPyg4AAAAAp0RlBwAAAHAATD1tfVR2AAAAADglKjsAAACAAzDx1I7VUdkBAAAA4JSo7AAAAAAOgGd2rI/KDgAAAACnRGUHAAAAcABGntmxOio7AAAAAJwSlR0AAADAAVDXsT4qOwAAAACcEpUdAAAAwAHwzI71UdkBAAAA4JRIdgAAAAA4JYaxAQAAAA6ARUWtj8oOAAAAAKdEZQcAAABwACYmKLA6KjsAAAAAnBKVHQAAAMAB8MyO9ZHsXK3kOHtHAAAAAKAKJDsAAACAA+CZHevjmR0AAAAATonKDgAAAOAAeGbH+qjsAAAAAHBKVHYAAAAAB2A08cyOtVHZAQAAAOCUqOwAAAAADoC6jvVR2QEAAADglKjsAAAAAA7ASG3H6qjsAAAAAHBKJDsAAAAAnBLD2AAAAAAHYGIYm9VR2QEAAADglKjsAAAAAA7AaO8AnBDJDgAAAACLFBUV6ejRozp06JAOHDiggwcPKjk5WSZT8RC8KVOmKCoqyqJzPfPMM1q9enWNrr9t2zY1bdrU4vYkOwAAAIADcPSpp9esWaOZM2eqqKjI3qFYjGQHAAAAQLWysrJqLdGJjIxUaGhote28vLxqdF6SHQAAAMABXC+zsQUHB+uWW25Rhw4d1KFDB7366qvat2/fNZ1zyJAh+vOf/2ylCP+HZAcAAABAte655x4NGjRIQUFBZfa7ubnZKaLqkewAAAAADsDRZ2OrycQAjoJ1dgAAAAA4JSo7AAAAgAMomb4Z1kOyAwAAAMCuUlNT9dVXX+nUqVMqKiqSv7+/WrVqpa5duyowMPCqz0uyAwAAADgAR19npzYtXLiwwv0Gg0F9+/bV1KlTFR4eXuPzkuwAAAAAdUxOTo5ycnLK7ffz85Ofn58dIqqYyWTS5s2btX37ds2fP18DBw6sUX+SHQAAAMAB2HI2tmXLllVYTZkyZYqioqJsEoPBYNAtt9yivn37qmvXrgoJCZGPj4/y8/N1/PhxffPNN1qxYoWKiopUWFiop556So0aNVKPHj0svgbJDgAAAFDHjB07VkOHDi2335ZVnaeeekoNGzYst79Bgwbq0qWLunTpovvvv1+PPPKIzp07p0uXLmnWrFn65ptv5OpqWRrD1NMAAABAHePn56cWLVqU+7FlslNRovN7ERERev31182vT5w4oW+//dbia5DsAAAAAA7AZMP/ric9evRQ165dza+3b99ucV+SHQAAAAAOrfRzOomJiRb345kdAAAAwAHU5amnqxMUFGTePn/+vMX9qOwAAAAAcGiFhYXmbS8vL4v7UdkBAAAAHIDJRGWnMkeOHDFvl67yVIfKDgAAAACHlZ2dre+++878ukuXLhb3pbIDAAAAOABbLipqb/n5+fL29q62ndFo1KxZs5SbmytJcnNz04ABAyy+DskOAAAAAJuaN2+eCgoKNHz4cHXv3l0uLi7l2iQmJmrevHnasWOHed/IkSPVsmVLi69DsgMAAAA4gOth/Zt777233L6UlBTz9hdffKEtW7aUOd64cWN98MEHZfYZjUZt2LBBGzZskLe3t2666SYFBwfLy8tLBQUFio+P19GjR8s8x3Tbbbdp+vTpNYqXZAcAAACARY4ePVrl8YyMDGVkZJTZVzIErTL5+fnavXt3pcfd3d310EMPKSoqSq6uNUtfSHYAAAAAB1CX1tl59tlnNWDAAO3du1eHDh3S6dOnlZWVpZycHLm5ucnf31/t2rVT165ddd999ykgIOCqrkOyAwAAAMAix44ds8p5/Pz81Lt3b/Xu3dsq56sMyQ4AAADgAFhnx/pYZwcAAACAU6KyAwAAADiAuvTMjq1Q2QEAAADglKjsAAAAAA7gelhn53pDZQcAAACAUyLZAQAAAOCUGMYGAAAAOAAjU09bHZUdAAAAAE6Jyg4AAADgAKjrWB+VHQAAAABOicoOAAAA4ABYVNT6qOwAAAAAcEpUdgAAAAAHQGXH+qjsAAAAAHBKVHYAAAAAB2BinR2ro7IDAAAAwClR2QEAAAAcAM/sWB+VHQAAAABOicoOrGLAe1uUnlNY4bEAbw9991i/MvvyL17WRzHH9V1cutKyC+Xh6qI/BPtrXPcwdW8VZIuQ4aSaNw/W7OenqX+/OxUQ0FDp6We0dt0mzX3hdWVlZds7PDgZ7jfYEveb8zNR2bE6p0x2srKydOjQIeXk5KhBgwYKDQ1Vs2bN7B2W0/PxcNWozqHl9nu5l73Nci4UadxnPynxXJ7CAn017NZWKiy6oq3xpzXpy5/1/D0dNPSWG2wVNpxIaGgrbd+2Vk2aBGntuo06dixeXbt01NTHH1b//nfqjt73KTPzvL3DhJPgfoMtcb8BV+e6SHYOHjyolStX6vDhw8rPz1fz5s3Vv39/3X///apX738j8TIzM/Xiiy9qw4YNunLlSplz3HzzzZo2bZoiIyNtHX6d4evhpkd7tau23aKf4pR4Lk9/bNtULw3pLNf//g6j8m/UyH/9qJe2HFLPkCA18fWs7ZDhZBa+NV9NmgRp6hPP6Z13l5r3v/ry83riiYma+8+n9diUZ+wYIZwJ9xtsifutbmA2NuszmBz8X/Wtt97Se++9V+Evv3Pnzvroo4/k7u6uc+fOaeTIkUpJSan0RjEYDJo2bZomTJhwzXEVLpl2zedwJgPe2yJJ2jC5T7Vt71n0rU7nXtDKh+5UeKBvmWOf7k7Uq9/H6tHb2mnSbW1rI9Trju+jn9s7hOtCaGgrxR2NVlJSitre2LPM3wEfH2+dStkng8Gg4Oa3qKCg4iGXgKW432BL3G/Wdbko1d4hVKpL8O02u9bu9O02u5Y9OfQEBZ999pneffddGY1GmUymcj979uzRc889J0maPXu2Tpw4Yf4D4OnpqaCgILm5uUkqTnRMJpNeffVVff/993Z7T87s0hWj/hN7Sh/uPK7Pdidq14kMXTGWTzwz8i9Kklo08Cp3rGTfLyfO1m6wcDp39u4pSfp2y4/lvvDIy8tXdPQueXt7qUf3zvYID06G+w22xP0GXD2HTXbOnz+v1157zfw6LCxMTzzxhObMmaMxY8bIy8tLJpNJX3/9tX788Udt2bJFBoNBHTt21Oeff659+/Zp+/bt2rdvn95//32Fh4dLKi4PvvDCCzIajfZ6a04rI/+i/vGffVq4/ahe+T5Wj3y5U0M++F67UzLKtGvo6S5JSs0uKHeOU//dl3w+v/YDhlNp1zZMknT8eGKFx4/HJ0mS2rQp/1wZUFPcb7Al7re6wyiTzX7qCod9Zmft2rUqKCiQwWBQr1699O6775qrNJI0fPhwPfDAA8rPz9ezzz4rk8mkrl27asmSJXJ3dze3c3V1Ve/evdWpUyf99a9/VXx8vNLT07V161bdfffd9nhrTune9i3VqUWAwgJ95eXuotSsAn2xN1n/PnBCU1bGaNmDvdSucQNJUq+wJlp9MEWLfjqmlwZ3lks9gyQps+CiPttd/Ic858Ilu70XXJ/8GhQPiczOzqnweE5OriTJ39/PZjHBeXG/wZa434Cr57CVnZ07d0qS6tWrpzlz5pRJdCQpPDxc48aNk8lkUkZGhgwGg2bPnl0m0SnN19dXzz77rPn11q1bay32umjybe3UrVWgArw95OnmqvAgPz3X/xaN7hqqC5eNeu+nOHPbv/Vqp6a+9bXlWLpGfLxNL3/3q+ZsPKD7P9oqv/rFv+f/5j8AAAB1RkWPbdTWT13hsMnOsWPHZDAYFBERUem00X36FD8MbzAY1LZtW4WFhVV5zp49eyooqHgNl8OHD1s3YFRo2K2tJUl7T54z7wvyqa/PxtyhER1bK7/osr7al6wdCb+p/43N9Mq9XSRJjbw87BEurmM52cXfbDZoUPE3m35+xd+MZmVV/M0oUBPcb7Al7jfg6jnsMLbs7OLFsUJDKx9/2qpVK/N2mzZtLDpvRESEtm7dqvT09GsLEBZp6FVcaSu8VHYq8ABvD83o214z+rYvs/+XE8XP90Q09bdNgHAax+ISJFU+Zr1NeIikyse8AzXB/QZb4n6rO+rSszS24rCVnaKiIknFs6pVpvQxPz/LxqkGBgZKknJzc68hOljqUFrxAmct/MvPvFaR9bEnJUkDbm5eazHBOW3dFi1J59KsAwAAIABJREFU6tvnDhkMZcdB+vh4q2fPrsrPL9DPMXvsER6cDPcbbIn7Dbh6DpvseHkVfzguLLRsvvjSi4ta0q4ujVWsbYnnclVYdLnc/tTsAr347a+SpIE3tzDvN5pMKqig/dexJ/X1r6fUoXlD3dWmae0FDKeUmHhCmzdvVUjIDfrbo+PKHJs9a5p8fLz16Wf/Zg0KWAX3G2yJ+63uMNnwv7rCYYexNW7cWHl5eUpLS6u2bU0Sl5KKjo+Pz1XHhrI2HUnTJ7sS1KllgIL9POXt7qpTWQXanvibLl42qldoY43t9r/nqS5cuqK739msHq0C1dLfWwaDQftTM3Uw7bxCA3z0ypAuqmdghgLU3JTHn9X2bWv15hsv6O67e+no0ePq1rWT7rrrNh2LS9DMWS/ZO0Q4Ee432BL3G3B1HDbZadmypRISEpSYWPX406VLl0qSmja1rBKQnJwsSQoODr6m+PA/XW8I0InMPB09k639qZm6cOmKfD3cdGvzRvpTRAv9KaJFmbK7m0s93XNjM+07lamf//uMzg0NvTXl9hs1qkuIPN0c9raEg0tMPKHukQM1+/lp6t/vTg24526lp5/Rm299qLkvvK6srGx7hwgnwv0GW+J+qxuMjDyyOoPJQcdzLViwQO+//74MBoO+/fZbtWjRovpO1SgsLFTXrl115coVDR48WC+//PLVn2vJtGuOB7CE76Of2zsEAACcxuWiVHuHUKk/NOlhs2v9+tvPNruWPTnsMzvt2/9vlq6SNXeu1XfffafLl4ufFenQoYNVzgkAAABYA8/sWJ/Djhfq3LmzevfuLYPBoLy8PKuc85NPPjFv33bbbVY5JwAAAADH5LDJTsOGDfX+++9b7XzZ2dkaOHCgBg4cKHd3d7Vu3dpq5wYAAACuFc/sWJ/DJjvW1qBBA40dO9beYQAAAACwkTqT7AAAAACOrC49S2MrDjtBAQAAAABcC5IdAAAAAE6JYWwAAACAA2CCAuujsgMAAADAKVHZAQAAABwAExRYH5UdAAAAAE6Jyg4AAADgAHhmx/qo7AAAAABwSlR2AAAAAAfAMzvWR2UHAAAAgFOisgMAAAA4AJPJaO8QnA6VHQAAAABOicoOAAAA4ACMPLNjdVR2AAAAADglKjsAAACAAzCxzo7VUdkBAAAA4JSo7AAAAAAOgGd2rI/KDgAAAACnRLIDAAAAwCkxjA0AAABwAExQYH1UdgAAAAA4JSo7AAAAgAMwUtmxOio7AAAAAJwSlR0AAADAAZiYetrqqOwAAAAAcEpUdgAAAAAHwGxs1kdlBwAAAIBTorIDAAAAOAAjz+xYHZUdAAAAAE6Jyg4AAADgAHhmx/qo7AAAAABwSlR2AAAAAAdgpLJjdVR2AAAAADglKjsAAACAA+CZHeujsgMAAADAKZHsAAAAAHBKDGMDAAAAHACLiloflR0AAAAATonKDgAAAOAAmKDA+qjsAAAAAHBKVHYAAAAAB8CiotZHZQcAAACAU6KyAwAAADgAE7OxWR2VHQAAAABOicoOAAAA4AB4Zsf6SHYAAAAA1EheXp5Wr16tTZs2KTk5WdnZ2QoICFBYWJgGDhyowYMHy93d3d5hkuwAAAAAjuB6WWdn9+7dmj59utLS0srsT09PV3p6unbs2KFly5ZpwYIFCgsLs1OUxXhmBwAAAIBFYmNj9cgjj5gTHTc3N91xxx0aNmyYunfvrnr1itOLY8eOaezYsTp9+rQ9w6WyAwAAADgCR5+NraioSFFRUSooKJAkRURE6N1331XTpk3NbeLj4zV58mSdPHlSZ8+e1bRp0/Tpp5/aK2QqOwAAAACq9+WXXyo1NVWS5O/vrw8++KBMoiNJ4eHheu+998zP6+zatUvbtm2zeawlSHYAAAAAB2AymWz2czWWL19u3p4wYYICAgIqbBceHq6hQ4eaX3/++edXdT1rINkBAAAAUKWkpCQlJiaaX5dOZipS+nh0dLTy8/NrLbaqkOwAAAAAqFJMTIx5OyQkREFBQVW2b9++vby8vCRJFy9e1P79+2s1vsqQ7AAAAAAOwJGHsSUkJJi3IyIiqm3v6uqqtm3bVtjflkh2AAAAAFQpOTnZvB0cHGxRn9KTFyQlJVk7JIsw9TQAAADgAGw58XROTo5ycnLK7ffz85Ofn1+5/VlZWebtwMBAi65Reqhbdnb2VUR57Uh2rpLnhFftHQLqiMvcawAA1AmXi1Jtdq23335bCxcuLLd/ypQpioqKKre/ZG0dSfLw8LDoGqXble5vSyQ7AAAAQB0zduzYCmdUq6iqIxVPMlDCzc3NomuUrLUjSRcuXKhhhNZBsgMAAADUMZUNV6tM6SrNpUuXLOpTVFRk3q5fv77lwVkRExQAAAAAqFLJNNJS2SpPVUq3K93flkh2AAAAAFSpQYMG5u2MjAyL+pRuV7q/LZHsAAAAAKhSSEiIeTs9Pd2iPqXble5vSyQ7AAAAAKoUFhZm3j58+HC17S9fvqy4uLgK+9sSyQ4AAACAKnXv3t28nZSUVO1QttjYWPN00x4eHrr11ltrNb7KkOwAAAAAqFJISIhCQ0MlSSaTSWvWrKmy/erVq83bkZGR8vb2rtX4KkOyAwAAAKBaI0eONG8vWbJEmZmZFbZLSEjQqlWrKuxnayQ7AAAAAKo1YsQINW/eXJKUmZmpiRMn6rfffivTJiEhQZMnTzZPO921a1f17t3b5rGWMJhMJpPdrg6nlpeXp9WrV2vTpk1KTk5Wdna2AgICFBYWpoEDB2rw4MFlVtYFrlZRUZGOHj2qQ4cO6cCBAzp48KCSk5NV8udtypQpioqKsnOUcBbx8fHasWOH9u7dq7i4OJ05c0ZFRUXy8fFRy5Yt1aVLFw0bNsxuD+PCuWRmZmrv3r06ePCg4uLilJKSojNnzqigoEDu7u5q0KCB2rVrp549e2rIkCFq1KiRvUOGk4uNjdWoUaNUWFgoSXJ3d1dkZKSCgoKUmpqqmJgYGY1GSVJQUJBWrFih4OBgu8VLsoNasXv3bk2fPl1paWmVtmnXrp0WLFjABwJckzVr1mjmzJllVmn+PZIdWMNPP/2k+fPnKz4+3qL2f/nLX/Tss8/abSE9OIdJkyZp69atFrX18vLSk08+qQcffLB2g0KdZ8nnvLZt22rBggUKDw+3YWTludr16nBKsbGxeuSRR8wzcLi5uSkyMlKNGzfWyZMntWvXLhmNRh07dkxjx47VypUr1bRpUztHjetVVlZWlYkOYC2xsbFlEp169erpxhtvVKtWreTn56ezZ89q9+7dysnJkSStWLFCiYmJ+uijj1S/fn17hQ0nUjI6Ijg4WJ6eniosLFRiYqJiY2NlNBpVUFCguXPnKiMjQ0888YS9w4UT69Kli9avX69Vq1Zp48aNSklJUVZWlho1aqTw8HANHDhQQ4YMcYgRPCQ7sKqioiJFRUWZE52IiAi9++67ZZKZ+Ph4TZ48WSdPntTZs2c1bdo0ffrpp/YKGU4iODhYt9xyizp06KAOHTro1Vdf1b59++wdFpzQzTffrOHDh2vgwIHlVgQvKCjQwoULtWTJEknSnj179Oabb+rpp5+2R6hwAt26ddPdd9+tnj17qmXLlhW2OXXqlJ5//nnt2LFDkvTee+/prrvuUocOHWwZKuoYHx8fjRkzRmPGjLF3KFViGBus6pNPPtELL7wgSfL399c333yjgICAcu3i4+M1dOhQ8zfyixcvtuvDa7h+nT59Wi4uLgoKCiqzf/To0frll18kMYwN1rF582bVq1dPffr0qbbt/PnztWzZMknF1e2dO3fK19e3tkNEHVZUVKT77rtPCQkJkoqHUZb8/xioy5iNDVa1fPly8/aECRMqTHQkKTw8XEOHDjW//vzzz2s9Njinpk2blkt0gNrQr18/ixIdSYqKipKbm5sk6dKlS9q5c2dthgbI3d1dgwcPNr+2ZIV7oC4g2YHVJCUlKTEx0fy6dDJTkdLHo6OjlZ+fX2uxAYAt+fr6qk2bNubXqampdowGdUXDhg3N2/w/FShGsgOriYmJMW+HhIRU+217+/btzbMUXbx4Ufv376/V+ADAlgwGg3n7ypUrdowEdUXJEDZJatGihR0jARwHyQ6spvQf2YiIiGrbu7q6qm3bthX2B4DrWVFRkZKTk82v7bnGBOqG+Ph4rVy50vy6b9++dowGcBzMxgaruZr/sZeepS0pKcnaIQGAXWzatMk8jMhgMKhHjx52jgjOqKCgQCkpKdqyZYuWLl1qngm1U6dOGjZsmJ2jAxwDyQ6sJisry7wdGBhoUZ/SQ92ys7OtHhMA2NqFCxe0YMEC8+v+/ftXOlkLUBPR0dEaP358pccNBoMGDRqkefPmydWVj3iARLIDKyr5RkmSPDw8LOpTul3p/gBwvZo3b555QoL69euzuCNsIjAwUPPnz2cZB+B3SHZgNRcvXjRvl0y5Wp3SK+teuHDB6jEBgC2tWLFCX331lfn1U089pZCQEDtGBGfStGlTjRo1SpJkNBqVm5uruLg4xcXFKSMjQ5MmTVKfPn00e/Zsi0dYAM6OZAdWU7pKc+nSJYv6lCwqKhV/AwoA16sdO3Zozpw55teDBw82fzAFrCE0NFSzZs0qt//48eOaPXu2du/erW+//VbHjh3TF198wfBJQMzGBisqmUZaKlvlqUrpdqX7A8D15MCBA4qKijJ/0dOrVy/Nnz/fzlGhrmjTpo2WLl2qjh07SpJSUlK4/4D/ItmB1TRo0MC8nZGRYVGf0u1K9weA60VcXJwmTpxofu7w1ltv1dtvv11mmC5Q29zd3fXUU0+ZX2/YsEGZmZl2jAhwDCQ7sJrS49LT09Mt6lO6HePaAVxvTpw4oYceesg8G2Xbtm21ePFiKtWwi44dO8rT01NS8UK2hw4dsnNEgP2R7MBqwsLCzNuHDx+utv3ly5cVFxdXYX8AcHTp6ekaP368zp49K0lq3bq1li5dSpUadmMwGOTj42N+nZOTY8doAMdAsgOr6d69u3k7KSmp2qFssbGx5mEfHh4euvXWW2s1PgCwlnPnzmncuHHmKaaDg4O1dOlSZsCCXV26dKnMmnck3gDJDqwoJCREoaGhkiSTyaQ1a9ZU2X716tXm7cjISHl7e9dqfABgDdnZ2Ro/frySk5MlFa9vsnTpUjVr1sy+gaHOi46OLjMbasn/k4G6jGQHVjVy5Ejz9pIlSyp9ODIhIUGrVq2qsB8AOKqCggJNnDhRx44dk1T8zfmSJUt45hC1onSVpjq5ubl65ZVXzK9vvvlmtWjRojbCAq4rJDuwqhEjRqh58+aSpMzMTE2cOFG//fZbmTYJCQmaPHmyedrprl27suIzAIdXVFSkv/3tb9q/f7+k4unyP/jgA9144412jgzOas2aNRo+fLjWr19vHvb9eyaTSdHR0XrggQd0/Phx8/4nn3zSVmECDs1gMplM9g4CziU2NlajRo1SYWGhpOLpMCMjIxUUFKTU1FTFxMTIaDRKkoKCgrRixQoFBwfbM2Rc5+69995y+1JSUswfDgIDA8s9S9G4cWN98MEHNokPzuHll1/WkiVLzK9vueUWtW/f3qK+HTp0qPA+Bary8ccf68UXX5Qkubm5KTw8XK1bt5afn5+uXLmic+fO6ddffzVPklHiySef1MSJE+0RMuBwXO0dAJxPRESEPvzwQ02fPl1paWkqKirStm3byrVr27atFixYQKKDa3b06NEqj2dkZJSbMCM3N7c2Q4IT+v2w3IMHD+rgwYMW9S0oKCDZQY2VXqvp0qVLOnLkiI4cOVJp+6ZNm2rmzJnq06ePLcIDrgskO6gVXbp00fr167Vq1Spt3LhRKSkpysrKUqNGjRQeHq6BAwdqyJAhLLoHAEAlRo4cqR49eig6Olr79+9XQkKC0tLSlJ+fL0ny9fVVcHCwIiIi1Lt3b915551ydeWjHVAaw9gAAAAAOCUmKAAAAADglEh2AAAAADglkh0AAAAATolkBwAAAIBTItkBAAAA4JRIdgAAAAA4JZIdAAAAAE6JZAcAAACAUyLZAQAAAOCUSHYAAAAAOCWSHQBOKSYmRu3atVO7du3sHYpDGj16tNq1a6e3337b3qHUWMnvNSYmxqrnvfvuu9WuXTutWrXKque1xNtvv6127dpp9OjRNr82ADgzV3sHAABVuXLlijZt2qStW7fqwIEDOnfunC5cuCBfX1+1bt1aXbp00eDBg9W2bVt7h1qlI0eOaMuWLfL19dW4cePsHQ4AAHUCyQ4Ah7V//349/fTTSk5ONu9zc3OTt7e3srKytHfvXu3du1eLFy9Wv3799Nprr8nd3d1+AVfhyJEjWrhwoZo3b+4QyU5wcLBCQkLUsGFDe4cCAECtIdkB4JC+//57TZ06VUVFRfL399eECRPUr18/tW7dWlJxxefw4cPavHmzli9frs2bN+vChQsOm+w4mpdfftneIQAAUOtIdgA4nOTkZE2fPl1FRUUKDw/XkiVL1LRp0zJtXFxc1L59e7Vv314TJkzQs88+a6doAQCAoyLZAeBw3njjDeXl5cnDw0MLFy4sl+j8nr+/v959912ZTCaLzv/2229r4cKF6tatmz755JMK28TExGjMmDGSpGPHjpU7fuDAAf3rX//Svn37dPbsWbm4uKhhw4Zq3ry5IiMjdf/995vjLj1JQmpqarlJE6ZMmaKoqKgy+zIzM7Vs2TJt27ZNJ0+eVFFRkRo3bqzu3btr/PjxatOmTbUxHz58WEuWLNGuXbt07tw5derUyfx+R48erV9++aXCa999991KTU3Viy++qD/96U/617/+pXXr1iklJUUuLi6KiIjQww8/rDvuuKPSf+OCggJ9+OGH2rBhg1JTU+Xt7a0//OEPeuihhxQZGVnmGn/+858rPU9NxcXFadOmTdq1a5fS0tJ05swZubq66oYbblDv3r01duxYNWrUqNrz5OXl6f3339fmzZuVnp4uT09Pde7cWZMmTVKHDh2q7Ltnzx59/vnn2rNnjzIyMuTu7q6QkBD169dPo0aNkre3t7XeLgCgGiQ7ABxKRkaGNm3aJEkaPHiwQkJCLO5rMBhqK6wyVq9erRkzZpiTK3d3d7m4uCgtLU1paWnatWuXgoODzR/iAwMDdeHCBeXl5alevXrlPmx7eXmVeR0dHa2pU6cqJydHUvFzSm5ubjp16pROnTqldevW6YUXXtB9991XaYybNm3Sk08+qUuXLsnHx0cuLi41fp8FBQV68MEHdeDAAXMMeXl5iomJ0S+//KIXXnhBw4YNK9fv3LlzGjNmjOLj483xX758WT/++KO2b9+u559/vsaxWGry5MlKTU2VJHl4eMjT01PZ2dk6cuSIjhw5otWrV+vjjz9WaGhopefIycnRsGHDlJSUJDc3N3l4eCgrK0vfffedfvjhB82dO7fC9200GjV//vwyCbSXl5cKCwt16NAhHTp0SKtWrdKSJUvUvHlz6795AEA5JDsAHEpMTIyMRqMkqW/fvnaOprzCwkLNnTtXJpNJQ4YMUVRUlG644QZJxclBYmKivvnmGwUEBJj7/PTTT1q1apVmzJih4OBgff/995We/9ixY3r00Ud14cIFDR8+XOPGjVPr1q3NydQHH3yg5cuX6x//+IfCwsLUvn37Cs/zzDPPqGfPnnr66acVFhYmSWUmerDEW2+9pfr16+udd95R79695ebmpsTERM2YMUP79+/XvHnz1L9/f/n6+pbp9/TTTys+Pl7169fXzJkzNWTIELm7uys9PV0vv/yy5s2bJ1fX2vnfT9euXRUVFaXu3burWbNmkqSioiLt2bNHr7/+ug4ePKhp06ZVOb30woULVa9ePb3xxhvq27evXF1dlZCQoNmzZ+uXX37R888/r5tuukkRERFl+r311lv65JNPFBAQoMcee0yDBg2Sv7+/Ll26pL179+r//u//dPjwYUVFRWnlypWqV4/VHwCgtvGXFoBDOX78uHn7pptusmMkFTt+/Ljy8/Pl5eWlF1980ZzoSMXf4v/hD3/QU089pd69e1/V+efPn68LFy5o0qRJmjt3rsLCwsxVmWbNmun555/X6NGjdfnyZS1atKjS84SHh2vRokXmREeSeXIHSxUWFmrp0qXq06eP3NzcJEmhoaFatGiRPDw8VFBQoB9++KFMn927d2v79u2SpH/+858aNmyYedKI4OBgvfbaa+rUqZMKCwtrFIulXnrpJQ0dOtSc6EjFlbfIyEh9/PHHCgwMVGxsrHbv3l3pOXJzc/Xmm29qwIAB5qQsLCxMH3zwgVq3bq3Lly/rzTffLNPn1KlTWrx4serXr6+PPvpIo0aNkr+/v6Tiylb37t31ySefqGnTpoqNja0y4QUAWA/JDgCHkpWVZd4u+bDoSEqqGJcuXSoTqzWcOnVKP//8s1xdXfXQQw9V2q5k+NrOnTt15cqVCttMmDDhqoaulda/f/8yyVKJRo0a6dZbb5VU/nmmjRs3SpKaN2+uIUOGlOtbr149Pfroo9cU19Xy9vZW165dJUl79+6ttF2nTp0UGRlZbn/9+vU1YcIESdL27duVm5trPrZ69WpduXJFt99+u2688cYKz+vj46M+ffqY+wMAah/D2ACgBm644QaFhoYqMTFRw4cP1wMPPKDbb79dbdu2vebkouQDuNFo1KBBgyptV5LgFBQUKCsrq8yQuRKdOnW6plgkVfkgfuPGjSVJ2dnZZfYfPnxYUvFwssqeoercubNcXV11+fLla46xIj/88IPWrl2rQ4cO6dy5cxVWkU6fPl1p/x49elR7zGg0KjY21vy65Hf3008/6bbbbqu0f0FBgSQpLS2t+jcCALhmJDsAHErpak5WVpaaNGlix2jKc3Fx0YIFC/TYY4/p1KlTeu211/Taa6/J09NTHTt2VN++fTV06FB5enrW+NxnzpyRVPxBOiMjw6I+lQ0HqygBqqmqZg0rGd71+4QlMzNT0v+SoYq4u7urYcOGOnv27DXHWJrRaNT06dP19ddfl4mzQYMG5mF4ubm5unjxYpXD6Kq650ofK3mv0v9+dwUFBeaEpioXLlyotg0A4NqR7ABwKKWnVD5y5IjDJTuSdOONN2rDhg3aunWrduzYoX379un48eOKjo5WdHS0Fi9erPfff7/cFNPVKZmYITAwUD/99NM1xXitVaZrZauZ8UpbuXKlvv76a7m4uGjy5Mm699571bJlyzITAUyfPl3r1q2zeJpyS5VU2x555BFNmzbNqucGAFw9ntkB4FC6d+9u/nD67bff1so1ShKBixcvVtqm9PMYFXF3d1e/fv30z3/+U+vXr9fOnTs1Z84c+fv7Kz09Xc8880yN4woMDJQknT9/3qLqgCMqmVa7pNJRkaKiIp0/f97q1/7Pf/4jSRo2bJgef/xxtWrVqtyMZ5ZUzH777TeLjpWeQjwoKEgSw9MAwNGQ7ABwKIGBgerXr58k6euvv1ZSUpLFfS39tr5BgwaSpPT09ErbHDx40OLrSlLDhg31wAMPmL/VP3z4cJkP9CUfuquKseQ5mytXrujHH3+s0fUdxc033yxJ+uWXXypts2fPnlp5XqfkOZySGH4vPz9fBw4cqPY8MTEx1R6rV69emeuU/O6io6OrTKIBALZFsgPA4TzxxBPy8vLShQsXFBUVVeU37VLxQ/JRUVHVVmNKlAwvO3PmTIUffs+dO6evvvqqwr5FRUVVntvDw8O8Xbqq4OPjI0nmhUIr0rp1a3Xr1k2StGDBgmrfj7Vng7OGe+65R5KUmpqq9evXlztuMpn0/vvv18q1S/6Njx49WuHxd999V/n5+dWeZ8+ePRUmPBcvXtRHH30kSerVq5f8/PzMx+6//365urrq/Pnzeuutt6o8f1FRkUVxAACuHckOAIcTEhKiV155RW5ubjp+/LjuvfdeLV68WCdOnDC3uXLlig4fPqw333xTffr00ebNmy0+f6dOncwr2D/99NM6dOiQTCaTjEajYmJiNHr06EorMP/5z3/0wAMP6IsvvtDJkyfLxLN9+3a99tprkqSOHTuaK0jS/55FysvL0zfffFNpbDNnzpSXl5eSk5M1fPhwbdmypUyl4LffftOaNWs0duxYvfrqqxa/Z1vp0qWLeTaymTNnatWqVeYE8fTp05o2bZp27959VRM4VOf222+XJK1YsUJffvml+bpnz57V/Pnz9eGHH1o0nbmvr68ef/xxbdy40VyBSkhI0MSJE5WYmCgXFxc9/vjjZfrccMMN5im1P/zwQz311FOKi4szH798+bKOHDmihQsXql+/fjpy5IhV3jMAoGpMUADAIfXp00fLli3TjBkzdOLECfOsZ25ubvL29lZOTo75gX6DwaA//elPFn+ArlevnubMmaNHH31USUlJGjZsmDw9PWU0GnXx4kW1bt1as2bN0v/7f/+vXF+TyaR9+/Zp3759koqf3fHy8ioTT+PGjTVv3rwy/Vq1aqXIyEjt3LlTf//73/Xcc8+ZP3iPGTNG48aNkyS1bdtWH374oaZOnarExEQ99thjcnFxka+vry5cuFBmFq+WLVvW7B/VRl566SWNGTNGiYmJmjFjhmbNmiVPT0/l5OSY/+0XLVqkwsLCMpWwa/XQQw9p06ZNSkxM1KxZszR79mz5+PgoNzdXJpNJI0aMUFFRkVavXl3leaZMmaIvvvhCU6dOlbu7uzw8PMxVNoPBoNmzZ6t9+/bl+j322GO6cuWKFi1apLVr12rt2rWqX7++6tevr9zc3DJrItljAgcAqItIdgA4rM6dO2vDhg3auHGjfvjhBx08eFDnzp1Tfn6+GjRooNDQUHXt2lX33nuvQkNDa3Tu22+/XZ999pkWLVqkvXv3qrCwUM2aNVO/fv00adIkxcbGVtjv7rvv1ksvvaSYmBgdPnxYZ8+eVXZ2try9vRUSEqK77rpLDz74YJkhTiXeeustvfPOO9q6davcuhVHAAABxElEQVTS09OVmpoqqfxkCJ07d9bGjRv11Vdf6fvvv9fx48eVm5srDw8PhYWFKSIiQnfccYf++Mc/1ug920pQUJBWrlypDz/8UBs2bFBqaqpcXFzUu3dvPfzww+rWrZtefvllSf9bpNUa/Pz89MUXX+idd97Rli1bdObMGbm4uKhbt24aMWKEBg0aZNHEEX5+flq5cqXef/99bd68Wenp6fL391fHjh01adIkdezYscJ+BoNBU6dO1YABA/T5558rJiZG6enpysvLk5+fn1q3bq1OnTqpb9++lZ4DAGBdBpO1598EAKAKycnJ6t+/vyRp69atCg4OtnNEAABnxTM7AACbKpmgIDw8nEQHAFCrSHYAAFaVkJCgf/zjH9q1a5fy8vLK7J8xY4ZWrVolqXgBTgAAahPD2AAAVnXkyBHdd9995te+vr66fPmyCgsLzftGjx6t5557zh7hAQDqEJIdAIBV5eXl6auvvlJ0dLSSkpKUmZmpy5cvKyAgQLfeeqtGjBihyMhIe4cJAKgDSHYAAAAAOCWe2QEAAADglEh2AAAAADglkh0AAAAATolkBwAAAIBTItkBAAAA4JRIdgAAAAA4pf8PE2Y78kfr1jIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x720 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set(font_scale=3)\n",
    "confusion_matrix = sklearn.metrics.confusion_matrix(y, y_pred)\n",
    "\n",
    "plt.figure(figsize=(14, 10))\n",
    "sns.heatmap(confusion_matrix, annot=True, fmt=\"d\", annot_kws={\"size\": 20});\n",
    "plt.title(\"Confusion matrix\", fontsize=30)\n",
    "plt.ylabel('True label', fontsize=25)\n",
    "plt.xlabel('Clustering label', fontsize=25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vis_data(x_train_encoded, y_train, vis_dim, n_predict, n_train, build_anim):\n",
    "    cmap = plt.get_cmap('rainbow', 10)\n",
    "\n",
    "    # 3-dim vis: show one view, then compile animated .gif of many angled views\n",
    "    if vis_dim == 3:\n",
    "        # Simple static figure\n",
    "        fig = plt.figure()\n",
    "        ax = plt.axes(projection='3d')\n",
    "        p = ax.scatter3D(x_train_encoded[:,0], x_train_encoded[:,1], x_train_encoded[:,2], \n",
    "                c=y_train[:n_predict], cmap=cmap, edgecolor='black')\n",
    "        fig.colorbar(p, drawedges=True)\n",
    "        plt.show()\n",
    "\n",
    "        # Build animation from many static figures\n",
    "        if build_anim:\n",
    "            angles = np.linspace(180, 360, 20)\n",
    "            i = 0\n",
    "            for angle in angles:\n",
    "                fig = plt.figure()\n",
    "                ax = plt.axes(projection='3d')\n",
    "                ax.view_init(10, angle)\n",
    "                p = ax.scatter3D(x_train_encoded[:,0], x_train_encoded[:,1], x_train_encoded[:,2], \n",
    "                        c=y_train[:n_predict], cmap=cmap, edgecolor='black')\n",
    "                fig.colorbar(p, drawedges=True)\n",
    "                outfile = 'anim/3dplot_step_' + chr(i + 97) + '.png'\n",
    "                plt.savefig(outfile, dpi=96)\n",
    "                i += 1\n",
    "            call(['convert', '-delay', '50', 'anim/3dplot*', 'anim/3dplot_anim_' + str(n_train) + '.gif'])\n",
    "\n",
    "    # 2-dim vis: plot and colorbar.\n",
    "    elif vis_dim == 2:\n",
    "        plt.scatter(x_train_encoded[:,0], x_train_encoded[:,1], \n",
    "                c=y_train[:n_predict], edgecolor='black', cmap=cmap)\n",
    "        plt.colorbar(drawedges=True)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Encode a number of MNIST digits, then perform t-SNE dim-reduction.\n",
    "x_train_predict = encoder.predict(x)\n",
    "\n",
    "#print \"Performing t-SNE dimensionality reduction...\"\n",
    "x_train_encoded = TSNE(n_components=2).fit_transform(x_train_predict)\n",
    "#np.save('%sx_%sdim_tnse_%s.npy' % (266, 2, 266), x_train_encoded)\n",
    "#x_train_encoded = np.load(str(n_predict) + 'x_' + str(vis_dim) + 'dim_tnse_' + str(n_train) + '.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEfCAYAAADsnan6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4FMUbwPHvXsmlkkZCFUJdSqjSiyiEjlSlShGRIqAoiPwUBUERC0oTkCJdio2iFAGpSrHQy1JCr6GG9Fz5/bEBjtACJrmU9/M897C3Nzv3zhrzZnZnZxSHw4EQQgjhSgZXByCEEEJIMhJCCOFykoyEEEK4nCQjIYQQLifJSAghhMtJMhJCCOFykoyEEEK4nMnVAQghhMj4VFXNCdQCqgBlgCJAXsAbiAeuAnuAtcBcTdMuP079ijz0KoQQ4lFUVf0FaJrC4lHA/zRNm5jS+qVnJIQQ4nFdAg4Cp4BowAsoATyNfvvHG5igqmouTdPeT0mF0jMSQgjxSKqqDgIigbWapoU/oEwIMAVomLTLAVTXNG37o+qXZCSEECLVqKrqBuwCSibtmq5p2quPOk5G0wkhhEg1mqYlAPOddlVIyXGSjIQQQqQ255F0Pik5QJKREEKI1FbSaftESg6QZCSEECLVqKpaCnjFaddPKTlOhnYLIYT4T1RV9QKKAi2AgehDuwH+AGakpA5JRkIIIR6LqqphwJqHFHEAC4AemqZZU1KnJKO0YQEqA+cBm4tjEUJkbEYgD/AX+rQ6TyoAyJGSgjNmzOCzzz6730fXNU27/h9iALgIdNc0bcXjHCTPGaWNWsBmVwchhMhUagNbnvDYAK5evUJAQIoKx8XFxT3zzDPuN27cSP7Rh5qmDX/U8aqqqkD/pLcGwA8IRZ+zDvSe0RKgj6ZpF1MSkySjtFEEOHrtWjR2e8Y+v4GB3ly5EuXqMLIsOb9pL7OfY4NBwd/fC/R7LseesJoQ4Di1asGZMw8vmT8/bNnC1KlTa48ZMyZ54f/UM1JVtTQwGT2xgt6eGpqmXXrUsZKM0kYIcPzKlagMn4yCgnyIiLjp6jCyLDm/aS+zn2ODQSEw0BugECkcBn0fIcBxQkLg5MmHlyxYEE6c+K/f90CqqlqA34EaSbu+0zSt06OOk6HdQgghUo2mafHA20672iUtP/FQkoyEEEKktq1ATNK2EX1A10NJMhJCCJGqNE1zAM6jI/wfdYwM7Xax8+fPMX3qdPbvOsRThfPQpFkTatasjZubm6tDE0KIJ6KqqhkIdNp19VHHSDJyocOHNZo3bkqRuFZcSrzJlj828NPcX8EjjpZtWtCkaVMqVarKli2biIq6SZ06z5E7dx5Xhy2EEI8SBjj/RX3oUQdIMnKhD98bQeWowZxx7MCfwnTkF3Yap7DJPpzDxzbyz8drOH7sEsGUxd9QiCHWIbRo1YKb12K4duU6T9coR4uWrShdOhRFUVzdHCFEFqWqaoCmaY/s3SSV9QU+d9q1U9O0E486TpKRC/2xbQNdHZ+wkY95i9NcYDf/+HzKsl96kjuPLwBrVh3gg4GreTF6BTmYwLKF06hAD/Z6LmSb9hdfz55Ovty5GTpoCKpaElUtIYlJCJHauqiq2h4YDyzVNC06eQFVVRWgHjAOKOX00ZCUfIEkIxfy9vTlSvwhfMiDG17sc5vJS69WuJ2IAOo3KsXEMZuZHB5KAtH04h++tTyH7//6kaftC8Tt38/Zfv3oM2QgHh4e+Hl4Mm3CFIoWLYaXlzdGo9GFLRRCZCFV0RfNS1RVdT9wGLiOPlouF/A0+rRGzoZomvZbSiqXZORCnbt1ZvmkKUTHX+IyGlbjTXz93e8pF+jrj5kSnGMH/zKDBK9Ezo8YzsUxn4PNRt5Ro/Bp3BiASyM+pHWbxoCCxc2N6tVr0bv3G1SrVj2dWyeEyEKc58wzA+WTXg9yBuivadqSlH6BDO12obcGDaJMw2DsxgRmUxdjbCCLZ+/Bar0zt+qpE1fYt/8UzzKMaC6zw3MaeUaPouSBAwS+8gqelSqRo0kTFEXh5po1sG4F02a/xI49Q5j7fReOHP2b1i80pWT50oSHP+lMI0KI7EzTtMlACfT56OYDO4ErQCKQAEQA/wDTgJZAocdJRJDBe0ZJizQ1QJ94NBTIC7ijj18PBzYB32qadvAx680FvAy0Qp9Gwxd9hu1d6Cf6J03T7KnTigczm81Mnj6FU6dOsnjxAr6f/zOXTkbwQrNvaP9SJa5ExLFw1k7qJX6OAzt4uJFz4Ov41KunV2C3YylW7HZ9MdMnM/z9+pSr8BQAIYVyMnZSO15oOwdaPk+NerXo270nb7wxEB8fH7m3JIRIMU3TNEADJqZF/RmyZ6Sqav2ka5L7ga+ANoCKvpa6GcgJVAEGAQdUVZ2WtLhTSupuCRwAPkmqIxh9yYcQ9Iz+PfC7qqr5UrNND1OgQEEGDRrCpm2bKFu6MueOWZnxyTF2fVOStlHrKOfozu+8h8WYA4+KFW4f51mpEjfXrcORmAhA3KnTlC6T9666CxYKxJGYSM4ePbAUL86UJT9QvGxxylQuy3cL5qVXE4UQ4qEyZDJCvxHmPBrDjt4tXAxMBZah3zi7pQewSlVVj4dVqqpqE+AH9HU/AKKApcC36NNX3FIHWK2qaorWBkktFouFH5f/zOixH+PwiOMQy9nmOYqxxnxcVY4QYn2GmD/uhOlZtSrmfPk43q4dN3//HUuePGz7I/yuOvftOYvJxxuDjw8e5crh37EjxsBAlG5def/LT/nuu+/Ss4lCCHFfGfoyHXoCmgos0jTtmvMHST2hYdyZkK8WMBK9t3SPpIn6FqCP/ABYBXR0rldV1dro67XnBEoDE4CuqdWYlDCZTLRv34l27Tqyc+c/HD8eTqFCr/H+Ox9w4vBWrk6MwODvi0/DhiSEh0NEBKYj1zk3cBCm/Pn4+OPfUAwKNWoX4eC+8wz78Dd8X38TgOg//yTX0KH4tmyJ9eJF/EeOZPinn1K//vPp2UQhhLhHhlxCQlXV1oBN07SlKSj7FTAg6W0CEKxp2j0rRqmqOgZ4K+ntYaCcpmlx9ynnvJyuHSiradr+x2xCCKm8hITD4WDbtj/57bdVrN/2J8ePaAQE56Jr+06M/ewrLKZ8eH78KqY8eYia8CVR/+7CzS8Hnq/0xqtGDS5PmoQ1IoKC8+dz9s03cQ8Nxb9dO07UrMnp4+fv+b7du3fy3sDXOXn0CF6+fvR+cxDduvVIlbZkJ5l9eYPMILOf46y2hMSTypDJ6HEkPe17iTtTT7TRNO2nZGXM6KM9bj3A017TtEUPqXMV0DDp7URN0/o/qOwDhJCO6xn9+uty+vbuTbxix6dRI9zLhRKzciUxBw5gs9sx+vjg27IlQX37Ert7N6f79KHo2rXE7NyJ/8yZrPtl7V31/f33dl5s2pCBOGjscPA38C7gW6gwnt7eFC1QkN69+lKtWo37xiPuyOy/KDODzH6OJRnpMuo9oxRL6gU591xC7lPsOe4kokjg50dUO9tpu+UTB5dOmjZ9nt379vPx8BHUctipfugIIzp25cjeI7Rt9QL2mzeJXL6cY02bcvq118jzySfE/Psv14YNZ/SwYffUN6R/Hz522PnQ4aAa0A/9RtuNE+GgXMPgdpqu3dpQu/bTjBwxnL1796R3k4UQWUym7xkBqKr6D1Ax6e1gTdM+T/b5J9yZkmK1pmmNHlFfHuCc067imqYdeYyQQshAK73evHmTdevWoGkHWb7uN06HH6NoiVK8++Yg2rdvc89flUXy+HPAZiO/0z4H4GVQWPfPYHLk8CAh3kr3TnO5ubsoVy0Heal7B0qW1qciCgtrgK+vX7q2MaPK7H+1ZwaZ/RxLz0iX0QcwPFLSErfFnHadvk+xkk7b/z6qTk3TzquqegHI7XT84ySjDMXHx4eWLVsD8M477z2yvKe7O0eio+9KRmcBNzcjHh761VA3i4n+b9VhZJ/dVI/+H1Onv0nNZ0rgsDt4a2BfnipagsZhDenx8qvkypUrDVolhMhKMn0yQn8GySdp24G+9npyxZ22T6Ww3tPcSUbqk4WWObXp8jK9p3zNOoeD/Ohj6F8GWrYpj9l8Z647i7uZOFsUmyzvsnhJT4oUCwbg0IHzvNRhDnOPF2Zy9YqEBIVQs051Br8zhMDAR64+LIR4UgunQnzsw8tYHvoEjMtk6mSU9FzRKKddP2iaduk+RZ0XebqYwuovOG0HPLBUFjRs+MccPXSQYuvXkctg4KLdjslkZGjFp26XcTgczJ3+N55xIdRqk/N2IgIoUSoPdeqX5N9ChbAE5uD85WOs+e0cixbPZsCAwfzz73Y2rN9AjhzetG/fmbfffg9FUTh4cD++vn4UKlTYFc0WIvObPwYizj28TFBeqNYgfeJ5DJk6GaFPVV4waTsWGPqAct5O24/4s+G+5bwfWCoLUhSFeYt+5uLFC6xZsxofHx8++mwUw4f+wsYNRylVOjerlx8mMtyfgo6qeOXYeE8d3l5m4r5fQOcXStKzd01MJiP7956j+0uf8WL7p9mw/S2uXYlmzKfraNNmC0cPnsZbycXNxAiKqUWZNnsaMTExGAwKhQoVkamLhMjiMm0yUlW1B/Cq065BmqYdfkBx56mwE1L4Fc6z1D5RvzbppmSGFxTk88D9oaH67bhXXunCpk2bWLx4Ebv+0jhyOIJS5upcdt/Jnh//pWe/Wvj66qfp6pVoVizfh4+vJ3361r6dSEqXyUuNmoVYtfBvFs3eRolCOek+sB7/G/gzTaMWUpIW2LCyYu9r1KxcE3eTJ7HWm5gsCp27duLLL7/MlMuxP+j8itQj5zjzy5TJSFXVBsAkp13zNU2b9KDyQBzgmbSd0t9mFqftlPam7pJRRtM9zOOMRCpVqiLDh+uDFi9dusTq1SsA2LcvN+1afEvLF0Kx2Rws/H4PlmeeJe+1o3f1aH5Y8DdHfj/Cd4k2KgHrDl/i5d4LwNvCRvO75ExU8SYXh+xLecY+kj9MQylSyodqNQuxdv2P5AyaxdIlvxEaWibVz0NayewjvTKDzH6OnUbTZWuZLhmpqloVfcoec9Ku34DujzgsijvJKKW9HOdyUSkOMJsIDg6mc+duADgcXWm1vT1Llixm7YbfiY11oBw9zsHwY1yOuEnOIB8cDgfTxvzO0kQbTyfV0RJIcMCIAE9KVffg2yXlqRw3mELU5aDHTNp2C6VP/2dYu/og2O3s3n2Ojp1fZM/OQ65qthAijWSqZKSqaiiwArg1Q/c2oLWmaY+69HYVfXZu0FckTAnncila+z27UhSFatWq317A79ixI/zzz998NOJDOrSazuuD6uLtbeHy9ZjbieiWZ4HeEVEs+vh53D3MrJg/m1IJL3EwfjGt2rxG27CJFLwey7MxCZzwMBMZm8jzzevRrm0X2rXriNlsTh6OECITyjQzMKiqWhR9zrhbI9v2Ak3utxb7fWhO2wVS+JXO5bQHlhL3KFKkGG3bdmD33oO88cb7LPvxLLOnaxgd+n80Z9uAwvn9AXi+ZVmslivsYxEoMPnTtbS6GMm6mATeB36LTeR1wB51hqnTRlOjnErHpvWZNHEsUVGZ9zKNECKT9IxUVX0KWMud536OAPWTz+T9EAeBFknbFR5WMOn7cjt9163jxWNSFIXu3XvSvXtPAPr26k7rn39gMfp6xRuAvu5mBr9VF4BzZ28ABjwIwOp2mbW/a0y33X3P7R3gq/3nCXY3805sIvkvX2bu3t3UnzqZZ+o3xNfPn+fCGlC5clVMpkzx4y2EIBP0jFRVDUZPRLeGcJ8GwjRNS+nzQgDrnbarJ02c+jB1nLbPPOZUQOIBJk6ZQYPefannZsEItPc00++jZtStX4KLFyL57OPVVLz5Dp1YSWK8DYfNgTVZHVcBxQEbYxPpi/4XxvdxcagXznNp7iwSJ3xFxxaNKVWyIAsWzE33NgohnkyG/tNRVVV/9Etzt2ZQuIieiFI6i8ItG9AnSM2BPmHqrRVdH8R5DaNHLmMhUkZRFIaP+IRhH47iwIH9dOjwPCOHr2D6N1s4deoauWzl8SYXCzzr0rrpi9gSrvPh8qV8Y7ejoE+vMdCgEGpQKGS9syq8gj6CZTb6esivAnXiYhk9+gPCw8M5cf4cFjc3Gjxbl7Cwhnh6et4vPCGEC2XYZJS0eN4KoGzSrmtAg4c8S/RAmqYlqKo6nTvrGY1QVXWZpmnxycuqqvoccGsiVTsw5bGDFw+lKAqlS4eyZ89x/vhjM0ePHsVggHUrNmI3rGJEh8E0bfo8169fo92B/VQ6d5Ya0dH8aTFxLcAT+/VY7Fb7Xd3609yZZqMcUM5spOjzJZj53QwM1WoSt2oDq5YvJiHBRp06zzF+3Ddcu3aNU6dOUrJkKXLlyn1voEKIdJMhk1HS5KdLgWpJu6KAxpqm/Ze1Cj5BX548B1AC+ElV1U6apt1evlxV1ZrAIvQ/tgHmaZq27z98p3iEmjVrU7NmbQA6d375rs/8/QNYuWk769evZf36dZxZ9gNeHiZuXo1hNPo07Ab0lRLHAAudjo0BQgoF4jCZcd/3N7MWdaVo8WAuXYzk9dcWU6lkYSwmE2Xc3dkZH0/9Bo0pX60GHh4etGzZBm9vee5DiPSUIZMR+vLh9ZzeHwA6q6raOQXHbtM0bV7ynZqmXVZVtSOwDP13WBPgtKqqa9FvRZQAnFeLOwA87qJ6IpUZjUbCwhoSFtaQkSNHs3v3Ts6cOcPE0R8x9expgoD9MTF8BNRMOmYFcMZk4PLlKBS7jWFDG1C0uD6yPzhXDgrm9KKk3c7chAQsCQlcB+otX8KvRzUcJjOD3htMr66vMvzDkTINkRDpJKMmo+Bk76skvVLCG7gnGQFomvarqqovAtPQh4h7c//F8zYBHTVNi0zhd4p0YDAYqFDhaSpUeJpmzZpz8OABrl27yrqVvzDi22lstVm56mZil1GhQdNQpk3disNsvp2IABLirWzYeJTj3Jliww8YDzSPjCR4wwbOvzOEb6ZPIvz4Mb74YqwsgSFEOsjwo+lSW9KS5KXQV9L+C3058njgJHqvqS3wnKZpZ10WpHgkRVEoVao0NWvW5oOPPmXNn//g3bUHZ0JC8M4ZxJVLvjgMbniWKMHG9XduM8bFJ+KwO+6axh0gHxB3/TrnKlbA/eefqGlwsPm3X6lRtjhd2rVix47t6do+IbKbDNkz0jStG9AtDeu/iH4P6ZO0+g6RvkJCCvHZZ1/etW/Bou8YPGQgY3ZZSUiwUrtOMQ7sPYcZ/YZkK6eysw0GzA4HvW5G8wlgSLBxCaiDA9/163hl0was3r4E5XuK/n378uKL7dOvcUJkA9muZySyjw7tOnJgj8aQd4YxbfJ22reewajJ/5JYtRovAe8pCsuAVxWFzzw9scXGMpI7/1MEo9+8vAbstNmIjbrBpSoVef3dwbzap4eLWiVE1iTJSGRpPj456NOnH5s3/o3Fy5cEzBg8PDA2asSXZjM9ihZlfu7cJPj54a3cO6V7bvSVbvMA9S0W3EuWoMjqVSxf+QtHj8qz0EKkFklGIlsIDg5m26btDOzQmQI3Ion74w88mzVDadECs6cnhXxycNMBW5MdNxuon7R9RVEweHpiCgrCp359Fiy47zgZIcQTkGQksg0/P39e7zeADSvWsXnNRroVLk7z6HhmjP6SLRu2MnHGHJoYjQxHf9isLbAZeBN9Pqp/Ae/nngPAHnkTPz8/F7VEiKxHcTgy9uJvmVQIcDyrLa6XHezfv4/xn41iz1/bOXPlMtVNZm4kJqB5ehI0YwZeVaoQu2cPpzp0ZP+uAwQEJB+Xd7egIB927z7E3LmzsNls9Or1GoGBOdOpNdlDZv8ZdlpcrxBw4gmrCQGOs+03iH/EWqAWD6jW4L9+X6rLkKPphHCV0qVD+Wb2dwBcv36NTZs28PXECcTs38PVSd9weex44nbv5KvPxz4yEQG8+kpvvl0wF2POQIw+Poz9ejwlixTDduMaJy5dpFzRYgz8cBR164alddNEdnBhOURHPLyMVxDQIF3CeRySjIR4AD8/f5o3b0Xz5q04cuQw8+fPwd/fn5fnLyZHDt9HHr9y5S/MmD8bvzKhxB84ABER+ATlJOrQfhagr2Xym3aIHi+1wxjgS2RkJKqq8tFHX1C1avU0b58QGYncMxIiBYoVK87w4R/xxhsDU5SIAMZ+PhZPm40OO3dyMCaGfXFxtDx3jgCgOuCJPv3HeGsiiZevUMIIlkMHefH5hkz9ZlIatkaIjEeSkRBp5Pylc5Sy25hos1EA/QL9bPTLEaudytUGTHYHO6MT2B5v5U/gs+HvceTIY09QL0SmJclIiDRSuHABmiWtxXSLAjwDOE8FvxUo4/S+AvCSzcZztSpRp3Yl/vnn73SIVgjXkmQkRBpp82J71t1n/3r0//HswCqgL/py6s4KANUc8OyRI7zQpB6//ro8bYMVwsUkGQmRRtq0acuJoCCGAzfRZ3IYAlwBRqBfrhuEPktvJafjEoAFSWUn2R387nDwZp9X+O67ubz8ckdee+1Vdu78Jz2bIkSak2QkRBrx9PRk019/8XPhIvijz3W3EfBAXzgrHv1yXQf05YwnAbPQ12UqDDRMqqcKkDc+jnffe5Ojx/4i4upOmj8fRs0KJZgwbgxXrlxJ34YJkQYkGQmRhgoWLMjabTv5edkqOrbrSMjzLQlt2Zq/LBb6u7nRSVGYi4E89GYMFRhkMlAOWMzd/3PGAS1al2PclPYc2XmG1gaFQWfPcXT0R1QpU5RGjeowa9YM5CF2kVnJDAxpIwSZgUHw4PMbHn6M1atXkpiYwK6/9rDm91VYrQm4e5rJHR/P31Y7tyYb+hnorMCm/e/xfr/vqbbhCO87/VyNAb7M4Y7N0w13iw/Llq0jT5686dK+jCCz/wyn6gwMS/qn7KHXlhP+6/elOnnoVQgXKFy4CH369Lv93m63Y7fbuXr1Kp3aNKWAptEEuGA285fVikdOT9zcTGzccozZyf7A6Q0MiYpn4rgX+PK95YRVLkvdhk3oM/Adzp8/y7FjR/Hz86dw4SKUKhWKp6dn+jZWiBSQZCREBmAwGDAYDAQHB7Nm81+Ehx9l8+ZN+Pn5UWH/HmbOmsyfW47h4WbkRoIN59ntIgGzQeHD3gsZHW+lHLDyl6U0/WUZnpYiGOMjiXVcIpfBwGWLhSEfjKDbKz1d1FIh7k+SkRAZUOHCRSlcuCgADRs2YePGdQzos4giT/kz+NhlFtnsmNCHh7/rZsQThR/jrdyaRKicw4EZByPijtAdA5/hwGy3cTg2hvojPsDTJwfHjmlcirhAtaq1admyDRaLxVXNFUIGMAiR0bm7u7Nq1SYmTZpJUL7SbPfJQQGTiZZAAaOBfcWCseMg+Wx2FZP+/Qw75qTt4sB7sTEMffM1Ll/fSlTMPj4d+iZVSxVi4cL52O32dGuXEM4kGQmRCSiKQtOmzVm0aCk7D59izsp1KC1bc9lkwM3fgwS7g/PJjvkXyAW3E9EthYCn8vpw5dQVzq88wOgbsYy+GcWUAX3p36NLurRHiOQkGQmRCZUrV4GpU2exY8c+1GJ1CQ7ORTejkctJn2vABKPCBTcjyWe4mwNYLWYObj7KxthEOgGdgb/sdjb8upwOHVpRvnwxnn22MtOnT5HekkgXkoyEyMTy5MnLp59+yeYde8nfriMFFIXc6PPbNe1ViwHvNiTMw8xUYA160lluNuLmaaZtgg13p7o8gRcdDqJiw5k+ty2D36/OvPnjqFi+NIPfGsTx4+EuaKHILmQAgxBZgJubG6PHfs1zjZvRs8fLeHgr1A5TKVMuPwUL5+THb7dy6ex19oVfZuXGN1iz6iCnD12EOOtd9ZwyKjR5vgwFQgIpEBLItLkvUb/6BA59506jJQ1YumI5JUqUdFErhaupqloKfWW+WkAokBdwB24A4cAm4FtN0w4+bt3SMxIiC2nYsDHjJ04mJsbK6JGriI6Kp1rNwoya3I6cRYNo9kIFgnPloHGzUH5VFH53OnYlsEFR2Ln5GFPGrufC+Rv4+nqgFi1AMfvzVI56mzf7D2DQkIEMeuctdu/e6apminSmqmp9VVX3A/uBr4A2gAr4oN+WzIk+c9Ug4ICqqtNUVfV6nO+QnpEQWUyLFq2oU+dZunRtS90aX1KiVD5OHL9CfFwcq9a/AUBAoBeff9OB9q8tIkdUPAaTgXNWO6rDQYMV+9nrZqTtjD/5YmoHwk+epwb52eo+npMnD3Lw4D942e0snjmdqlWrM++HZTIsPOt7Gijl9N4O7AaOoM8BnBt9dZRbE4f0AEqoqtpA07TYlHyB9IyEyIL8/PxZtnQNW7bs5O2Bn/PL8t+pULESS3/efbtMleohlKpVjEsF8uJbLh9tjAp/2xz0ASYl2Jgek8iQnguw2uNYy/844bGdoLg4tickcMlq5TRg2b6Vd9/q98A4RJazE+gD5NQ0raKmae00TeulaVoLID/wuVPZWsDIlFYsc9OljRBkbjpBxjq/x4+H88ILTcid14MSpYLY8McpovzzYFRLYvpxMcujE+5aysIBBBoUOr7+LHN+OIhy8QrfJibSxqlMBFBAUWjVuRu9X+lFyZKlSG8Z6Rw/iVSdm+7kNrDGP7ykyQIFqz3W96mq2hqwaZq2NAVlvwIGJL1NAII1TbvxqOPkMp0Q2UShQoX5449drFixnDnzZnPx7A28zQHc/HkZHjGJJP9tEQ/YTQaatyrL5Imb8PUPoFTE3ZNwBqHfvf7p7+0sWbaEKeMm0bhx03RqkUjusLKNROXhv/fNii/FqfZY9Wqa9tNjFB8OvAa4Jb3qAY88Xi7TCZGNuLu707r1iyz56RcO7D7Ewq8m8s/Wf+nc9w3eNRiISSrnAD41KpQtk5eEBBsWvxy4VazIUkW5q75/AcXhoOOBAyjXr9GzZxd2796V3s0SGUhSL2i/066QlBwnyUiIbCpHDl/Kl69IcHAw7743jMLNWvCUyURrg0I5TzcW5PPj7U+aM/Kjtfi274D3W2/xkbs7HwF7ge+A1sCXwLfoiclktdGyRQNu3ox0YctEBuB8f8KYkgMkGQkhMBqNjJ8+m5/WbuJ81aoHAUHaAAAgAElEQVSEKwYUb3fatp7BiVyl8O/3OpaiRTGGhfGpYqG22czHwGSgW1IdKtDQYiRv/hyMGzfGZW0RrqWqqgUo5rTrdEqOk2QkhLitVKlQli39ja1bd9O3zzAsHj4QZ+Pq3LlEvNID09//4Obrh61oUVorCo2THW9TFIqXyMWBA3tu74uNjeXkyRPExcWlb2OEq7RBf/4I9B7S7w8pe5skIyHEPXLnzkO7dh3564+/eb1qTRpfusLQVi+yc9tOlv/6K4nnzjHZZOKi0zG7gQ1WO1evxFChQhXsdjsfjxhJaIkSNHqmKSWKFaF3z1e5ceO6q5olnEydOjW/qqohyV5+jz7ywVRV9QBGOe36QdO0Syk5VpKREOKB/Pz86flqbz77+DM6dHgJDw8PihUrzoihw4lWFIoCPQwKndxN1LGYqBGmsn/vRXr1eo1xX41l2YwtNImdQYw9EWNIcVadPEbp8iVZsHC+q5uW7S1cuHAzcDzZa8BDD3q0cUDBpO1YYGhKD5RkJIR4bC93fYW//z3AS7378qNPDhYm2IhMsHL6lIPfftuMj08Opk6eQljseJZaepNz0sc8tWIxBb+bS46uXRn2Zj9a1q3FwoXzsdlsrm5OttS+ffva6M8aOb/GPml9qqr2AF512jVI07Tkk8Y/kDz0mjZCkIdeBdnn/N76PaIkDf22Wq3kz5eTNo6FrCk/gbw/zsJhtxPRqRPq3r0Mjo0lAfjc05PC9RsxftqsJ/7uzH6OU/Oh18OnxpJofcRzRiZfihcY8F+/7y6qqjYAfuHO8lnzNU176XHqkIdehRD/mZLs+SOTyUTxQmU5Hf4nhgB/AKI2bSJg3z42xMbe/sXTJCaG4qtXMGLE+/j5BRAW1oBSpUqnc/Tiv1BVtSr6Q623EtFvQPfHrUcu0wkh0sQHH73HAfd5RP65HmtEBHFbt9I5Juauv4A9gbC4OH5a/C2TPx9Jo2er82yVsoSHH3VV2OIxqKoaCqwAbs3QvQ1orWlawuPWJT0jIUSaCAtryLzF83hjwJuEN22Ke9myHDSZwHr3GkobDAqlrkUzJtFOMDD9xAka1alOkzZt8PPLiZeXF4GBOWnWrAXBwcGuaYy4h6qqRdHXbAxI2rUXaKJpWvST1Cc9IyFEmqlWrQbbt23np5nzCcuZi2UOB6uSPnOgL4xz1e5gWaKdsujrEAwFmsfHs2n9L8ybN5VjJ1ezZes8qlUrQ+/er7Bt25/IvW7XUlX1KWAt+n8y0JeSqK9p2rUnrVN6RkKINFetWg2qVavBHx070+vVbnjExhJrTeSqzUpNNwMeMYl3lW8MLDsfyYK1/SlYKBCAQwer0rHNdFasXEKefE+xesV6/Pz8XdCa7E1V1WD0RHRrCPdpIEzTtIsPPurRpGckhEg3NWvWZvvew3y9dAWTf1qOw2Lh3wQbicnK7QDyKbBl05Hb+0qUzE2NWoXpN+A5ihX3oFmzutJDSmeqqvqjX5ornrTrInoiOvVf65ZkJIRIV0ajkbJly1O5cjXGjZtMrM1BF+AC+uI3M4F5QHMHnDtx9a5jTSYjgTm9+GJcG25GXWH79q3p34BsKmkZ8RVA2aRd14AGj/Ms0cNIMhJCuEzz5q1o/EI7VqPPrOkNzAF+BTYoUKBo0O2yx45GsGPrcWrXKcbJE1cpWCiAHTu2uSbwbCZp8tOlcHshpCigsaZpex581OORe0ZCCJfq1KkLq3/4kVpYGeZw4At8ARwwGdk/Zh0njl8hOjqedb8d4o236zHkrZ/QDl0gONiHL7/8hNjYGwwePPyeZ51EqhqJvkjeLQeAzqqqdk7Bsds0TZv3qEKSjIQQLlW9ek2Kl63E2r072Bbgh8PhwFylCmZV5caMGSya/w9ms0LHLlXYsfU4efL6MmpMK2aM28CVU9eYOnYc//y5g6mzv5MBDWkn+Zj6KkmvlPBGv/L6UHKZTgjhUoqi8OOyJdQNa0RkTAzUqYMjMJCYefN4ve8AypevQJnyT7Fv7znWr9UY8HY9+neajWPxv2yIjONvm52ntv5B5dBibNy43tXNEU9IekZCCJfz8PBg3rwFnD17huXLl5CYaKXRsjcpVqw4TRo0pFnzMLz8vHF3N7Nr5xmMZ68zK8HGrQtzs4FyDgdtO73InBlzaNiwiSubk+VomtaNO+sopglJRkKIDCNfvvz07t3vrn1lypSjeo26bNixFQ+TG2tWHSAsNhHnO0QK0DwxEa1wYXr17cmw94bTrFkLgoKCyE5M+VsB1oeXyaC/9mXW7rQRgszaLZDzm5omT57IhMkTibx2icp2B1sS71564hmziR1mIxXK5+P8+RtcvhRF1669GDZsZIYe3JCas3Yvsi4hiofPxuONF+1MLf/r96U6uWckhMgU+vTpx4E9h/hl+ToOGs0MA6KBGOAjReFfm43q1Qty/vwNrNdisEbHM2fSeEoUzsfVq1cfUbtwNUlGQohMpXz5Cqz98x++9vTET1HwMxr5MigntRqV4tiRCK6dvU7M9VheBUYDodFRVAstyrVrkpAyMklGQohMJ3/+p9i66yB5i6kkKgrRkde5djWaAE83vBNszAEmAH2BjUCY1UqViqVktF0GJslICJEp+fn58/eWHRw9dIg3+g/ksBbBlQs3cABNncopwADAzw169uzMuXNnXROweKiMOawiHaiqmgPoCryAPhNJAHAJOAgsBL7TNC3edREKIVKiSJEiDBr0P06cDGf5D4uxoI8nMzuViQQCAr2oVrUIixcvYMCAQa4JVjxQtuwZqapaC30hqPHAM0AewAI8BTQAvgV2qKpa0mVBCiFSTFEUvp44nbAmzbGjTyd0SzTwvtHA812rkP8pHyZOGkdwcA6KlirCN9OmyMzfGUS2S0aqqlYEVgEFknYlACuBGcB6wJ60vyywTlXV/OkepBDiiXz77TwGjxzNJ4qCCrRSoIDJQEiTUtRrUJI5EzcRd/0GJgXcom/w8agP6PZKFyIjb7g69GwvWyUjVVXdgB+5s177P0ARTdOaaJrWQ9O0ukAZIDzp8zzA/PSPVAjxpHr1eo3D565SICyMv/IF8uqQBjzXvAwv1ptAk6h4zgLRDvgkNgFzfALr1q0ktGJpfvrpe1eHnq1lq2QE9EJ/OAzgCvp67WecC2iadgB4Hrh1v+gZVVUbp1uEQoj/zGg0Mm/eD/xvyCj++jOWaV/vhah4ZgA5ATegB9ABsMcl4N2tGwPeeYvTp//zGnHiCWW3ZPSa0/bnmqZdul+hpIQ06wHHCSEyAYPBQNu2HVi4cDmv9XmL8oAxWZnqVjtF8/uRZ8dKDAYHs2bNcEWogmyUjFRVLQ6UcNo1+xGHOH8epqqqd+pHJYRID6VLl+Ev7lzuuGWtu4mmHSqxYH4nunerxMyZU6gYWoHhQz/g5s1IV4SabWWbZAQ857StaZp24RHl/4Lbkzy5A9XTJCohRJorXlylzNOVaQrsAc4AHyiwwdON1u0qAtCtR3US4m00uDSbP2dG0LBufdasWcWNG9ddGXq2kZ2SkfMw7X8fVVjTNCv68O/7HS+EyGTmL1lBwc7dqKUoFAP+rFmYmUt64uvnCYBBUVBQsGPlkHktxy6E81KPrpSsUIovvvrctcFnA9kpGRV32k7pXcrTTttqKsYihEhnFouF0WPGs23fUcw+nngF+WCz2Rna/3vCKoymefUvMCkKcz1a4NG6PkXXrqXob7/h06oVYyaOlamE0lh2SkaBTtsXU3iM86W8gFSMRQjhIkFBQdjtDg7sP8eLYRMpu3I/f0bGsfBGHGpCHBY3yD18OOa8eTHnzUvuYcMwh4TQb0B/YmNjXR1+lpWdkpHzAISU/kQ5l5MBDEJkETmDAikXmpd2Cnzo0Bf2qYn+1DtRUVgv3vl7VVEUvGrW5GrEDapXrswLL7Rk5sxpMnNDKstOc9O5O20npPAY58E3Ho/7hUkLZmV4QUE+rg4hS5Pzm/Ye9xy/1qcf4z8azgSr/a79OYBQm41Tv/9OQMeOt/dHbdyIYjBy2dvMjYtn+XvUB3z86SecPHYEPz+/VGiByE7JKM5p2y2Fx1icth+7fy4rvQo5v2nvSc5xly69WDBrNtsPHqSV0/5E4DAQ/fFIfMLCMLi7c3nKFBJPn8GnYRjWzevo2KkiT+VT+fnH3aglirJ92z68vLwe8E2P5rTSa7aWnZJRlNN2Sns5zuWiHlhKCJGpGI1GJs2YR4Nalalit9MSuAYMdDNSokQuDpy4ypHq1VGMJgLNpTFgIGHLemZ9257SZfIC0OqFCvTsNo/mzRqxZt1GDIbsdNcj9WWns+e8zGOuFB7jXE6WiRQiCylatBjBISH0MBnwNRsoaDYS1aAkY+d0oXP36pjdzZSytaFd3ELAQXAuH4oUC+LM6WvExyeiKAqdulThxNkj9OvXx9XNyfSyU89IAxolbRd4WEEnzuW01A1HCOFqNWvWYZ9/IhOndcTd3YTFoq+ClDevL54hBTl0eR37r/6AMQEiD93k2bKj8DYZiTcodO5Rg6B8vriVCuXnFcsYb/0ak8m1v1KLGGoSj/2hZSwZtA+SnZLRQaftCo8qrKqqCX0G7/sdL4TIAgYPfpdKleYSH5eIr69+Vd7hcPDjD7vw7NIT/yJFuNiuHZ0UhTCHgz+BhTYr84D3pv3JSS8L3qPGcP3g/zh16iSFCxdxaXu+t5/gyj2THt0tEAuVDEHpFFHKZadk5PzEmqqqai5N0x72vNHT3FlqIg7YmmaRCSFcIleuXPTu3Z8Xm0+l7xvP4h/oycKFuwhP9CZ3y5ZcffVVRgCDkoZxt0OfimU88E1cIvUNZjyqVIG4WAIC5FHE/yJj9tfSgKZph4FDSW8VoMsjDunqtL1O0zQZwCBEFvTeex9SuWpdvhizkQ+/2sHJ59qTe+5CDBYL13bs4JVk5bsBa4FSQEx8PBeGDqVW7Tr4+fmne+xZSbZJRkkmOW2/rapqzvsVUlW1BPCy066v0zQqIYRLzZoxh8ZhTYg8G4FnteoYPPRLdm4m0z3TtUSgPwH/C+Bht2PctYuZU75N54iznuyWjL4BTiRtBwErVFXN51xAVdWS6D9ntx6S3aRp2sp0i1AIke4MBgOTJ03jpRfbc7xlS463aMXRemHYrDYGmd1uP6SYALwFVAF6GRQSLCaWLfwZb295Tui/ylbJSNO0BKANEJO0qzJwTFXVX1VVna6q6lpgH3DrLuQF4KX0j1QI4QpffP4VB3cf4v2OnZnwv/c5ePgkbnXrkcdgoK7ZQG6DwgY3I7bqIYyf3xWjQeH48WOuDjtLyE4DGADQNO1fVVUbAfPQh25bgCb3KboXaK9p2un7fCaEyKICAgLp0aPX7fffzF1EteqhlGxaiDebhVK4yJ2RaPnzB7J40UIaNWrqilCzlGzVM7pF07TN6MO23wC2oPeAEoCzwBrgFaBy0vLjQohsrlBIMa5eibsrEV29Es3pM1fYuHkdJ0+ecF1wWUS26xndomlaJPoIzfGujkUIkbENHTqSxo3r4B/oTv2GJflkxEoO/XUKN7sDtwAz//vfW3z33U+uDjNTy5Y9IyGEeBylS4fy9tvvM3vGVjq3mU74jpPUtjtQgZxXo9mydi1HjhxxdZiZmiQjIYRIge7dX8VuA2ucFbekNZCGoz8D4g5069LOpfFldpKMhBAiBby8vChatDh+6AloItAYeBv4Cbhw7CjR0dGuDDFTk2QkhBAp9MYbg4gFOiTb/yz6L9Pnnnkm3WPKKiQZCSFECrVo0RoDkPx5j0j0CSxPnD7K1atX0j+wLECSkRBCPAZTTj/eVLg9TVA8MMCo4OttAZORdl3aEx//8Jmzxb0kGQkhRAodO3aE2JgEgqsVopjZyNNuRnIZFJYCEVHxeBuNHDys8eW4Ma4ONdORZCSEECk0d+ZcisS3ZP+Bq3h6mLlps+Nhd7Dc5sAGrI+PJ8/168yY8Y2rQ810JBkJIUQKnT19gUK2MLxuBPJGZBzeNgdzgJro69JUAr4H7NevuzTOzEiSkRBCpFDN56py2ON7LnCSN4DDQNVkZSoB0Q4HiYmJ6R9gJibJSAghUujFF9tjDT4J2EkEQoGNycpsAUKCgjGbzekeX2aWbeemE0KIx+Xl5cWKtatoUa8OI08dZyjQC3CgP2u0DX0J6YGDhrgkvhcNIcRjf2gZSwbtg0gyEkKIx+Dr68eiX1dTtazKakWhnN1Od/RnjdwVhZtGI4WLFHVJbP87E8s5q/WhZfKaTKwpkE4BPYaMmSKFECIDy5UrN1YPbw74+7OpcAjXTAaMgX4Y6z2HYrGwb99eV4eY6UjPSAghnoDDZMBii+f1zlVo2boDlyOiGPXJOv416Qv0ZUWqqroB5dBXXq+a9G9x9MGEAB9qmjb8SeqWnpEQQjwBxZ5IixahdOpSFS9vCwULBTJufCsUm5Vt2/5wdXipTlXVLuhXI3egzxPbGVC5k4j+E0lGQgjxBDxMBqpWKYDD4eDy5ShiouNxdzdTrFgwPj4+rg4vLQQAlrSqXC7TCSHEEygUUpSlP+3i67HriYiIwmq18Ww9lSPaBSZ//Yqrw0tLp4HtTq9Pger/tVJJRkIIkUIOh4MdO7YzafyX7Nu9G6PFyFcT21L72WJERsYxesRKFIOBaVOm8fqbr5MvX35Xh5yafgAWapp2wXmnqqoJqVG5XKYTQogUcDgcDB30Bv3btqDymlX0dzjwtNo5d/IqiqLg6+vBiE+aY1AUfli4imdr1eb06VOuDjvVaJp2JnkiSk2SjIQQIgW2b9/Gmh8WsTs2lg+BL4F/rXYmfLaWiEs3ATC7mchfwJ+4vBCt2Hi9fz+XxpyZSDISQogUWLV8Cd1iY8nhtC8EaGRQ2LzxCAA3rscQHn4Fn3r1wNONrf9ul7WNUkiSkRBCpEDE5Qhi7rP/htXG9WuxbP8znG6d5uDbujW5330X7zp1UNwtHD6spXusmZEkIyGESIFiJUvzDXcvOf43sMEGkyds5K0Pfieq02sEDh8JgFfNmihAQECAC6LNfCQZCSFECjRr1pwYi4USJhNtjEYaAs8oCjY3E3aDkcBJ0/Dv2AnFoP9ajd27l2D/gKw2oi7NyNBuIYRIgaJFi2F1ODD4ePOHBxATT+z1WPK6G8Fi5vyQwRT4dhYGHx8iV67kxoIFLPh5havDfqCpU6fmHzPmnuXRr2ua5pKVASUZCSFECgUEBXHj4gUMMQpz4q3UA3Zdj6WLJZFrkRpHnnkGR2IiPv4BLF28lAoVKro65AdauHDh5vvs/hAYns6hAHKZTgghUqx39154WG1MiLcShj4pWwVgUbwVYhN4as4cHMCSBT9QpUryNWAzlvbt29cGCiV7jXVVPNIzEkKIFOrf93W+GvnBPXPflEr691Snl2je4HlCQ8umd2iPrWfPnmd69ux5wtVx3CI9IyGESKFdu3Zi4t6lxneh95KmjPua6d/OTP/AsgDpGQkhRAr17fsKNi83+sVbsVjt1Af+RV9qvE79xrRq9YKLI8y8JBkJIUQKOBwOLlw4S40qBTm89zx9o+K4nmgnp4+F64k26jZo4OoQMzW5TCeEEI9w/Hg4gwa9gT06ntDNxxh5OYpn4qx4KfDc86EkWu2Z4j5RRiY9IyGEeIivvvqUKd9MIMjfjc6AyWrnf0AHoI3Vzry5f2GwmKhYsZKLI4W3AwKIcTgeWsZTSZWFWVOdJCMhhHiANWtW8/Wkr3jrnXr8Ons7JYEJwF7AP6nM60CleCujRn7Mex8MdVmsAAPX2jkV/fBkVMBLYW/LdAroMchlOiGEuI/+/XvTvXsHKlZ6irWrDnIk/DK/AK9wJxEBlAZqAl9PGUdk5A2XxJoVSM9ICCGSmTt3Fr+u+Illv/UlX3499Xw3ZzvjP1xJnfuUt3t44DCZ+fbbaQwYMCh9g01Hqqruus/uok7bvVVVTd7vOqdpWpNH1S09IyGESGb8pHG0e6ny7UQE0LFLVQIK+jMBuOxUdhew1eEg1/Bh/Lw6485Fl0rK3efl5fR5rvt8XooUkGQkhBDJREVF4uPtds/+osVzEalAYaCn2Uxbd3dqWizk/OILPMuXJ+LSpfQPNouQy3RCCJFMaGgFFn23g05dq5KYaGP5T7s58PcpNm04gtXdA6uHBz82b45boULkb9QIU86cXJk5k/Llyrs69DSlaVqaDcWTZCSEEMmM/mg0z9avSdNG32C7EkntRBthVjsOk4k18fEYatTg5tq15HrnHWyRkdxYvpyIMWN499e1rg4905JkJIQQyRQpUpQSxUtzTDvI23FWPkh6dqef1cpHwOiNG/Hr14+r8+eTcOoU9ugYurZ/idDQMq4NPBOTe0ZCCOEkOjqaAQN6c+zwAeKio+mf7CHS/kCc1UrEN1OJ2b4dIm/S7YV2fPrpPQvViccgPSMhhHAycFBf4hPCWb2xHw2rjyEqwXbXc0VRgBtgtDpYsmYjJUuWxmw2uyjarEN6RkIIkeTixYusW/sbwz5uzPlzN1ASbbwH2JM+twPvAM0Bkz2RsmXLSyJKJZKMhBAiycWL58md15/o6AT6tJtJTwecBkoCLyf9uxJ4Bsjh4+PKULMcSUZCCAGcP3+OLz75ghPh5xk1eAl5YhKoCPwOzARqJf1bFZiKQqeefVwZbpYj94yEENleRMQlqleqRrzBisFgYf3mY3wILAE6AjWSXlfRV3n18/fjtf5vujDirEd6RkKIbK9dmxeJN1oJ7NuHvHNng6LQEdgHvAr8AfwMVAPcn32WczExXLp00ZUhZzmSjIQQ2VpcXBwHTx0h6N3B5OzbG68KFQhs0IBxwGb0ydYGAIPRR9KZ8ufHp0o1Zs6c5sqwsxxJRkKIbC0mJhq7LRHPp5/myjvvcLZyZRJ372Yy0BJ9BJ074Al8D8T+8gu2K1cxGOQuR2qSsymEyNb8/QMwmN241KkT3W/epI/NxlVgIHAMfcrpt4DngXDAkZBAXPgRevX62YVRZz2SjIQQ2dqsmTMxRCdQ1xHPOKf964Cn0HtD3QEj8BGQEB/P633fIigo2AXRZl2SjIQQ2ZbVauWzUaMp6qhMa7bc9Zkn8CwQAPQATAYj5z09+X7+YqpXr5n+wWZxcs9ICJFtRURcIjHeyin+YEeyz+zAHqAJcALoM3k6Ow+flESURjJsz0hVVTegDhAGVAJKAIHoPyOXgd3AKmCOpmk3n6D+6ui979pAPsAGnAHWADM0TduXCs0QQmRg/v4BWOOv0gwHi9F/4bQBooFh6CPprEBQYBCtWrVxYaQpMym3N3HxDy/jbkmfWB5XhkxGqqp+jt4z9ntAkaeSXs2AkaqqvqZp2sIU1m0BxgK97/OxL1Aa6Keq6khgpKZpjvuUE0JkAQcO7MeKnbHovZ+O6NP+GICmwLvAy4qBz8aMdV2Qj2H6fDsXIx5eJlcQ1K6W8S6KZchkhD5wxTkR3QS2A2fRezAl0J8/MwD+wAJVVYM0TZuQgrpnAh2c3h8GdqCfi1pA/qTtDwEz8P5/aokQIkOJiopiz55d/LJ0ObPnzsKBPnS7OnAU+B8wHX32hTVuboz4YhxNmjzvwoizh4yajAASgR/Qfy42appmc/5QVdUSwDzg6aRdX6mquknTtN0PqlBV1Z7cSUQ24HVg8q3ej6qqJmA48F5SmaGqqq7XNO331GmSEMKVJk34ijGjP6YwDs4mWvEjgLic3nx5OYKR6CPmPgMswNcmE/tPXsRoNLo26Gwioyajhej3gsIfVEDTtEOqqoYBO4EQ9J+jd4F29yuvqqo78IHTrlGapk1KVqcVPQEVBF66VQ69FyaEyMRGjvyA7yaM5R+gGPpfox9zjc+uOPjKy4ut0dE0AlYDfwK+AXkkEaWjjHfhENA0bfjDEpFTuevAp067Gj2keAv0gQoA19CTzIMMQf9ZBaiqqurTDykrhMjgjh07wtwpExmGnohA/+t1KHY8HA5iEuxsCAhgqNf/27vz6KiqPIHj3ywkgQDNIiAQcGnkh06jgjY6iijiBj22Gy5n9KAScW+1223UccSFcWvQwXZDRERQUdwXFA+tIo0ti45DH+wfimwiCCYYFkkgSc0ft7BuQir1sr6q1O9zTo63Xu599TuVkt+7792lLXOzcuibcQklm7azbt13IUadXpIyGdXRAq/cXkQ6xal3mleeqaql8U6oqutwc952O6MB8RljQjZt2jPktcqgd7XjmbgRc/mDBtHzz3+m4w3Xk5PfiRzyycxoRVHRjyFEm56S9TZdXVQf7RavXz3UK38c4LwfASfV0NYYk2JemjmDdl3aMn3DVk4rK//l+ApgRUYG+zw2kay2bQFoc8QRLDnjbLIrIvTt2y+kiNNPS+gZ9ffKO3BzkKoQkQ7A3t6hzwOc169zYP1CM8aEbe3aNWzdtpUJj5/L0m7tOCc3mzeBR3Ej6NqcfvoviQggT4TMzh24uLCQvLy8sMJOOy0hGV3klf8aZ15Q32qv1wQ471qv3FFEutQ1MGNM+DZsWE9FRSV9DujKs29eRq9rjmX8gALeOrYPm1tlkXfUUVXqV5aVUVlSwqWX1TQV0TSVlE5GIjICONE79Ficqp298pbanhd5NlR7He9ZlDEmifXvfwitW+cw76Ovadsuj9GXH8Njsy7ht8OE1t32YuNDD7Fzrbv2jOzaxcb772fAgAEUFPQKOfL0krLPjESkK/Ckd2iuqr4bp3pbr7wj4FtUr9e2xlq16Ny5zk1C0aVLu7BDaNHs8216tX/G7TjuuBO56Y+vcu2fhtH/kB4sWLCSp6cspKwig5zuBaw8+Xe06tmDyuJNHD5wIK+/+or93ZpZSiaj6OTUmbjVEsBtTX9RLU38G787A75N9RWeWgds94uiom1UVib3akJdurRj06Y6L+1nArLPt+kF+Yy7dO9J7hnn8PQXP1L+2nyyDzyIHi+8zKbHH2frO7PJjeRw7P79uL8XI3sAAA8oSURBVGP6i/Tp4wZ/N9ffLTMzI2UuXJtSnZKRiFwNXN3IMQxX1ZV1bPMUbnV3cOsYXqCqtU0I8G/L5QR8j+rLCQbtURljkkynDh3J2VRK54f/UuX4rtWrGTpkCBMmPEKPHj3jtDbNoa49o70AaeQY6rSGbHQR1YuiLyPAGFWdnaDZNq8ctIdTvd62GmsZY5JSUVEREyY8wKyXp7P1p63QKovynFZ0u/k/iEQibJ09m9z165n65hwbNZcEUuo2nYjcCtzgHbpRVacGaFrslduLSF6AQQzdajmHMSZJRSIR7r13LE888SgZpTvpiLtPn7ergh8mPcWKWa/QpksXckvLeHn6S5aIkkSdkpGqjsUtJNrsorcIx3mH7lHV8QGba7XXvXGrddfGn6y9WVUTLMxujAnT4sULufW6q1iz4huIVJKbk0U73IKVl+GGDv8NOKm4mJ777sdb731Efn5+qDGbmJQY2i0iFwITvUOPqGrgrR2ia9j5Q7UHBGjm1/kq6HsZY5rfZ599ysgRJ3DqcuXjigqmVEboUFpOK+AKYv/QHY3bKG1VcRF333d3aPGaPSV9MhKRs4CngYzooWeBa+txqg+98rEB6vt1PoxbyxgTuluuuYIrgbuA3+AWovwI2IjbtdXXF8grKOCFF6YTiST3aNd0ktTJSEROAZ4ntt7cLKCwnruvvuGVz41uKRHvfXtQdTLt6/V4P2NMM/j5559Zt/LbKishA+wDtMcNvd09+qgSmJqRwc+7dlG6fTuVlZXNGaqpRdImIxE5BniV2FDs94Dzq2+yVwdv4HaKBbeaws211L2XWAL8TFUX1/M9jTFN7LxzzgdgqXdsJTAQl3yewe0dcy1wMvDNPvuwc+1afn3Qv9h+RUkkKZORiAwE3iY2vHoecKaqBp2wuofo6Lm7vEO3i8jlIrL79h8iki0idwGjvHq31vc9jTFNa9u2bcx7901uAO4BPsXN9zgN+Hfge+BLYBEwA/gqK4uMwYOp3LqVsTffFu+0JgQZyXjPVEQ2Av7CpNOAoNOh/0dVv67l3C8A53mHlgOf4UYWDgb8BanGqep/Bnxf377ASluBwdjn23RWr17F8ScPpXRzMe9GIhThbndsx91O+Y7Yg2aACcBTGRksb9eOAfsdyOwPPggh6j15KzDsB6yq52n2BVYumwU7E8yIzGkLB42koe/X6JJ1nlH1FbJH1VirZrOAuMkIN2G2BDfaE9zzzOqrepfjLrTuwhiTdNasWcXxRw6kfaSSYZEIFwOHAP8AXgMmUTURgVs77IdIhMMO6s8br77VzBE3jzk3QMnq2uv8ap9fklFSSdZk1GRUtQy4XESmAaOBY3C3lCtwF1MfAJNV9R/hRWmMqc1pJw5lZEU5k3HPGsqBc4A/4fYoWoy75N83Wj8CTAa2ZWRy8aiLyM5Ou3/6kl5S/kVUtfpFTVO8xwKqbllujEkBa9asZuPmIh4g9tA7Gzcj/re4FZQzgUHZ2fxXeTk9gSczM1nYtSsVFZUsWbKYkSPPDSd4E1dSJiNjjIln9ux3qGTPPV3a4xa63AhMBa7KzOSeYcPIKCmBk0+m+7nn8vOIEfTq1RuTfJJyNJ0xxsSTm5tL9w6teara8SeBf8PNySgEBu8qJ/P44+k0cyadRo9mx9KlVPxYRGHhpc0es0nMekbGmJQyfPjvuP32m7glA5ZE4EhgDvAFbg7IbgMjlcx74CHKi4up+HYdP81+g3N+fxa5uXXaKMA0E0tGxpiUEYlEeOD+e6msjJAZcZNa3wGW4Ca9do7WKwdeJoNWZWUUT/wLuVk5XHPpH7jlNptblKwsGRljUkIkEuHMM09h8eJF/KpDHltKdjBvZwXzcYufXgRcH607jkzK2rXnG/2W0tIdtGmTT2amPZVIZpaMjDEp4b777mbDD9/w5pyr6FnQkXXfbWb02U9zwMZtCG34mh3MJ0JWVg4HHz2YjydPJTs7m7Zt24UdugnAkpExJiVMnz6FRyadRc+CjgD0LOjIw5PP55JR09ix5VgqKv/GS3Pe5NBDB4YcacsmIu2BC4GRwAG4tT434rbaeRF4Pjqfs06s32qMSQnFxT9xgHStcqxP365sKSnlm8o5nHPBSEtETUxEBuMez00EhgDdcSPqewEnAVOAhSJyYF3PbcnIGJMSDj74N3z816obNM/7cDn5bXO4e9x/M37CwyFFlh6iC1i/R2wX7J3AbNx+cx/ixpMAHAzMFZGCupzfbtMZY1LCbbfdw5gx57OlpJQBh/fmiyVreHDc+1x77Q2MGXNF2OG1aCKSA7wC7N6nfQlwuqp+59U5CHgL2B/XY5pBsI1MAesZGWNSxJAhx/Hcc6+y4OMdXFn4Cm+/to6pU19m3LhxYYeWDi4jttRfETDCT0QAqroMOBXY/bxoiIgMD/oG1jMyxqSMQYOOYPr0V8MOIx1d6ZUfVNWNNVVS1WUiMpXYrghX4m7lJWQ9I2OMMXGJSF+gn3fo2QRN/N+fICLVlxGskSUjY4wxtRnqlVVVNySovwi3xyFAHm5Xj4QsGRljjKmNP0z780SVVbUcN/y7pvZxWTIyxhhTG38n7DUB26z1yhKkgQ1gaBpZ4Pa2TwWpEmeqss+36aXyZ+zFntXQc7UPMLNnd51JkyYVjB8/vvqvf1LVn6od6+yVfwgYin8rr1OQBpaMmkZ3gI4d8xPVSwqdOwd6vmjqyT7fptdCPuPuwIp6tt0CbB49n45BKpeWlpZOnjz5kxp+dScwttox/8PdETAev16gP44lo6axCDgGWA9UhByLMSa5ZeES0aIGnKMY6IPb8DahGTNmUFJSUtOvqveKwA1C2G1nwHj8telaB2lgyahplAHzww7CGJMy6tsj8hVHfxIqLCyksLAw6HlLvXJOwDb+DoaBelM2gMEYY0xttnnlQL2cavW2xa3lsWRkjDGmNn5vq1vANn69QL01S0bGGGNqo165d9xaVfn1NG4tjyUjY4wxtfnKKw9IVFlEsoH+cdrHZcnIGGNMbT70yiIiiW7VHUZsq4lS4NMgb2LJyBhjTFyquhz4Z/RlBjAqQZMLvfJcVbUBDMYYYxrFY175RhHZq6ZKItIPuNg79GjQN7BkZIwxJpEngVXRchfgXRHp6VcQkQOBt4lNkp2nqoH2MgLIiEQijRCnaWlEZBWwTx2arFDVPk0TTeoQkfa42xQjgQNw63JtxD3EfRF4XlXL4p/B1ERExgJ31LHZMapqk88biYgMBD4B2kQPlQFzcSvN7IvbamJ3B2cDMEhV1xKQ9YyMaSQiMhi3dP5EYAhuiZdcoBdwEjAFWBi9gjQmpajq58ApxFbuzgVGAIXAMGL5ZCkwrC6JCGw5IBPMNGBrgjqbmiOQZBW9anyP2Ciinbirxu+B/YFjcf+zHgzMFZFBqvpdGLG2AIuAhQHqfd/UgaQbVf1ERPoDFwFn49bD64T7/38Zrvc/oz69f0tGJog7VHVV2EEkKxHJAV4hloiWAKf7yUZEDgLewiWm7sAMXIIydfeuqo4NO4h0papbcL3/iY15XrtNZ0zDXYa7Zw5QBIyo3utR1WXAqcRWMx4iIsObLUJjkpwlI2Ma7kqv/KCqbqypUjQhTY3Tzpi0ZsnImAYQkb5AP+/Qswma+L8/QURaxK5wxjSUJSNjGmaoV1ZV3RC3prMI2B4t5wH/2iRRGZNibACDCeIwETkD6IHbubYI+F9ggapur7Vly+cP0/48UWVVLReRpcCRXvsPmiKwFqybiIzCzePKBzYDK3GTLNfU2tIkLUtGJohZcY5vF5GngLtUdXNzBpRE+nrloP8QriWWjKRxw0kLl0d/9iAic4HbVTXQ4pwmedhtOtMQ+cB1wOfRuQfpqLNX/iFgG/9WXqdGjMW4yZefiMh1YQdi6saSkYlnF25ezKXAIUB7IAc3R+b30d/tti9uraq9mznGZOAPQNgRsI1fzwYwBLcMuBP3nK4b7vvYARgE3AtsidbLAh4SkQvCCNLUj92mM/EcqapFNRzfgEtEb0Xv2z+Du6gpAO7DzcxOJ3leeWfANv7s9NaNGEtLNjHORNcS3KCQRSIyGZhN7NbpIyLytqr+1EwxmgawZJRCRORq4OpGPu1wVV1Z/WCcRFS9zjQR2Z/YApYXiMgtqrq+kWNMZqVeOSdgm1yvHLQ3ldZUtThAnW9F5FTc2mi7e01jgAebODzTCCwZpZa9aPwH3rmJq9TqQeAm3BV+FnAC8FxDg0oh/sZhQXs5fr1AG4+ZYFR1uYi8SGwDuFOwZJQS7JmRaZDo0O6/e4fSbUVq/4o90XbMNdVLeMVv6myuV06372PKsp5RConeMx8bchg18UeHdY5bq2VS3NU3QO+Abfx62rjhGNL7+5iyrGdkGkMbr5xuk2C/8soDElUWkWzAHwb/Vby6pt7S+fuYsiwZmcZwqFdOp8ELAB96ZRGRRLfqDiO21UQpYJMzG59/UZBu38eUZcnINIiIDKXq9uTzwoolDKq6HPhn9GUGsQfn8Vzoleeqqg1gaETRnuf53qG0+j6mMktGZg9BV5IWkU7AE96h5bg5H+nmMa98o4jsVVMlEekHXOwderRJo2oh6riy+X3Ar73XMxo5HNNELBmZmvxdRMaLyCHxKojICcBnVF2b7WZVrWzy6JLPk8CqaLkLbjWKnn4FETkQeJvYJNl5qjq72SJMbdeLyPsicpqI1DgVQUQKROQ54Hrv8OuqOr95QjQNlRGJRMKOwSQZEVlF7NbbBuBL3LprO3BrqQ2i6q05gLGqemczhZh0RGQg8Amxh+dluCHG63HLJQ0ldvG3ARikqmubOcyUJCJjiU2s3gH8H/Atbvmf1rgLosOpOjr4S2BIdItskwJsaLdJZO/oTzybgD+o6sxmiicpqernInIKMB03dDsXGFFD1aXAeZaI6q01cET0pyYVwBTgj7a9SWqxnpHZg4j0Ao4GjgIGAl1xqz+0w60YsBFYDLwPzFTVsjinSjsi0h63Pt/ZQB9cT3ITbpHPF4EZ9nnVjYjk47bcOCr6396472MnoBy3n9EyYD7wrKquDilU0wCWjIwxxoTOBjAYY4wJnSUjY4wxobNkZIwxJnSWjIwxxoTOkpExxpjQWTIyxhgTOktGxhhjQmfJyBhjTOgsGRljjAmdJSNjjDGhs2RkjDEmdJaMjDHGhM6SkTHGmND9P2isJMzugEpzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualize result.\n",
    "train_new = False\n",
    "n_train = 2660\n",
    "predict_new = False\n",
    "n_predict = 2660\n",
    "vis_dim = 2\n",
    "build_anim = False\n",
    "    \n",
    "vis_data(x_train_encoded, y, vis_dim, n_predict, n_train, build_anim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import keras\n",
    "from keras import backend as K\n",
    "from keras.models import Sequential, Model\n",
    "from keras.layers import Input, LSTM, RepeatVector\n",
    "from keras.layers.core import Flatten, Dense, Dropout, Lambda\n",
    "from keras.optimizers import SGD, RMSprop, Adam\n",
    "from keras import objectives\n",
    "\n",
    "\n",
    "def create_lstm_autoencoder(input_dim, timesteps, latent_dim):\n",
    "    \"\"\"\n",
    "    Creates an LSTM Autoencoder (VAE). Returns Autoencoder, Encoder, Generator. \n",
    "    (All code by fchollet - see reference.)\n",
    "    # Arguments\n",
    "        input_dim: int.\n",
    "        timesteps: int, input timestep dimension.\n",
    "        latent_dim: int, latent z-layer shape. \n",
    "    # References\n",
    "        - [Building Autoencoders in Keras](https://blog.keras.io/building-autoencoders-in-keras.html)\n",
    "    \"\"\"\n",
    "\n",
    "    inputs = Input(shape=(timesteps, input_dim,))\n",
    "    encoded = LSTM(latent_dim)(inputs)\n",
    "\n",
    "    decoded = RepeatVector(timesteps)(encoded)\n",
    "    decoded = LSTM(input_dim, return_sequences=True)(decoded)\n",
    "\n",
    "    sequence_autoencoder = Model(inputs, decoded)\n",
    "    encoder = Model(inputs, encoded)\n",
    "\n",
    "    autoencoder = Model(inputs, decoded)\n",
    "    autoencoder.compile(optimizer='adam', loss='mse')\n",
    "    \n",
    "    return autoencoder, encoder"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
